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Heat Transfer to High-Quality 
Steam-Water Mixtures Flowing in a 
Horizontal Rectangular Duct’ 


E. J.. DAVIS? and M. M. DAVID®* 


Heat transfer and pressure drop were investigated 
for steam-water mixtures flowing in an electrically 
heated horizontal duct of rectangular cross-section. 
The investigation was primarily concerned with the 
region of high vapor fractions and mass flow rates, 
where two-phase convention-controlled heat transfer 
predominated. 

The experimental data were found to be in good 
agreement with those of other studies in which 
circular tubes were used. The data in the convection- 
controlled region were correlated by equations devel- 
oped from two flow models, a separated-annular 
model and a homogeneous model. 

Two phase pressure drop data were correlated 
satisfactorily by the Lockhart-Martinelli correlation. 


I the vaporization of a liquid stream flowing through a closed 
channel, three regimes of heat transfer can occur—nucleate 
boiling, forced convection-controlled heat transfer, and a regime 
of liquid deficiency or burnout, which occurs beyond the point 
of maximum heat transfer coefficient. The first regime is 
characterized by bubble growth and nucleation at the heat 
transfer surface, and the heat transfer coefficient is a function 
of the heat flux. In the convection-controlled regime the heat 
transfer coefficient is independent of the heat flux, as in single- 
phase convection heat transfer. The third regime is characterized 
by a liquid deficient condition at the wall and occurs when the 
vapor mass fraction becomes so large that the liquid film on the 
wall is removed. 

Nucleate boiling has been the object of many pool boiling 
and flow boiling studies in recent years, and numerous theoretical 
and experimental investigations in this field have been reported. 
Two-phase convection-controlled heat transfer has been studied 
much less extensively; and only a very small amount of work 
has been reported in the liquid- -deficient range other than ge 
boiling burnout measurements. Much of the early work 1 
two-phase flow heat transfer was concerned with situations 
where both convection and nucleate boiling contributed to the 
heat transfer process:**:4.5©, and the work was generally 
uncoordinated with other investigations and limited to specialized 
studies and to the measurement of average overall or integral 
coefficients which shed little light on the mechanisms involved. 
More recent two-phase convection studies have measured and 
correlated point coefficients and avoided other complications of 
cennnconccnseusscccccucneccesccccescessoncesneccencescsnccecsencnccacnacencaces 
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earlier studies. Dengler“:®, Silvestri and coworkers, Fikry, 

Kvamme™, Rohsenow and Clark“2%, and Mumm*®) have 

studied the steam-water system; Johnson and Abou-Sabe®, 
Groothuis and Hendal“”, and Fried“® have investigated the 
air-water system; and McNelly (19.20) and Guerrieri and Talty@ 
have studied various organic two-phase systems. However, the 
published correlations and data are still limited in breadth of 
applicability and are often inconsistent. 

The present investigation, although extending into the 
nucleate boiling and liquid-deficient ranges of operation, was 
particularly concerned with two- -phase convection-controlled 
heat transfer with the steam-water system. The study used 
high quality (usually above 30% steam) mixtures, flowing in 
a horizontal duct with rectangular cross section, which was 
heated on one vertical face only. 


Apparatus 

The overall apparatus is shown schematically in Figure 1. 
Metered water and steam streams were mixed by spraying 
water into the steam to produce a mixture of the desired quality 
and at the desired pressure. This mixture was passed through 
the heat transfer test section, and then led to condensing tees 
where the mixture was condensed and cooled. 

The enthalpy of all entering and leaving liquid streams was 
measured by calibrated chromel-alumel thermocouples installed 
in the flowing streams, and the steam enthalpy was measured 
by a throttling calorimeter. Flow rates were measured by 
calibrated orifices. These measurements plus that of the heat 
input to the test section permitted mass and heat balances to be 
made around the system for each run, and the : enthalpy balance 
was usually found to be accurate to within 5%. 


The heat transfer test section, shown in Figure 2, was 
designed to measure heat transfer coefficients in essentially a 
differential length, which has a rectangular cross section 0.7 769-in. 
high by 0.260-in. wide. The test section consisted of an 8.25-in. 
long converging inlet portion, a 7.25-in. calming portion, a 6.0-in. 
heating portion : heated electrically on one vertical face, a 2.50-in. 
unheated portion, a 3.00-in. window section, another 3.50-in. 
section, and a 4.00-in. expansion section. The window section, 
which had the same internal dimensions as the rectangular 
heated section, was installed to permit visual or photographic 
observation of the flow patterns. Pressure taps were located 
on the wall opposite the heater area and were connected to seal 
pots or manifolds which led, through 1-in. needle valves, to 
calibrated pressure guages and a manometer system. The 
pressure measurements permitted determination of both the 
temperature of the two-phase mixture in the heated section and 
the pressure drop along the test section. 

The rectangular duct was constructed from 26-gauge type 
302 stainless steel whose thickness was accurately measured 
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Figure 1—Schematic diagram of the apparatus. 


before the duct was fabricated. A thin-walled duct was used to 
minimize longitudinal heat transfer from the heated area, and 
heavy steel back-up plates were clamped on, but insulated from, 
the vertical sides to eliminate the possibility of bowing at the 
higher pressures used. To permit obtaining the inside wall 
temperature of the thin-walled duct, provide an accurately known 
heat transfer area, and ensure a uniform distribution of the heat 
flux, a 0.100-in. layer of pure copper was electroplated on one 
vertical face of the 6.00-in. heated section. Five calibrated 
chromel-alumel thermocouples were mounted in 0.032-in, diam- 
eter holes drilled vertically along the center of the copper plate, 
at the positions indicated in Figure 2, and the inside wall 
temperature was determined by correcting for the temperature 
drop across half the thickness of the copper and across the duct 
wall. The high thermal conductivity of the copper and the 
small size of the thermocouple holes make the accuracy of this 
method essentially dependent upon the accuracy of the thermal 
conductivity data for the stainless steel, for approximately 90% 

of the calculated temperature drop from the thermocouples to 
the inside of the wall occurs in the stainless steel. Thermal 
conductivity data from a number of sources were compared, 
and the data published by Smith®® were used in the present 
work. The thermocouples were calibrated before and after 
installation and are estimated to be accurate to within +1.0°F. 


The heat flux was supplied electrically by means of a stainless 
steel heater strip with the same area as the copper plate, with 
heavy copper leads fastened to the two ends. The copper leads 
were connected to an 18 K.V.A. transformer which supplied 
the power at low voltages and high currents (less than 10 volts 
and up to 400 amperes). The heater strip was clamped to the 
outer surface of the copper and separated from it by a layer of 
isomica 0,0009-in. thick. The heat input was determined by the 
power supplied to the heater strip, corrected for losses along 
the leads and to the surroundings. The test section was insulated 
with glass wool and the remainder of the installation with 
standard pipe insulation. 


The heat losses from the test section could not be determined 
by an overall enthalpy balance because the heat input was only 
a small fraction of the total enthalpy of the flowing stream. 
The heater was therefore calibrated, both by heating a stream 
of water and measuring the temperature rise and by measuring 
the power input with stagnant air inside the test section, with 
the test section at the operating temperature. 

The extent of temperature clevation (above the temperature 
of the two-phase mixture) of the duct walls adjacent to the 
copper-plated heating section was investigated both by rough 
calculations and experimentally by thermocouples mounted on 
the outer surface of the duct, as close to the copper strip as 
possible without actually touching it. Both methods indicated 
that longitudinal conduction of heat (horizontal conduction, in 
the case of the top and bottom of the duct) through the duct 
walls, away from the heated section, was negligible. From the 
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Figure 2—Heat transfer test section. 


calibration studies, it is estimated that the heat flux was known 
to +3% 


Experimental Data* 

Experimental data were taken over a range of total mass 
flow rates from 50,000 to 600,000 Ibs./hr.-ft.*, qualities from 
30% to 90% steam, pressures from 25 to 150 p.s.i.a. ,and heat 
fluxes from 60,000 to 260,000 B.t.u./hr.-ft.2. The upper flow 
rate limit for the 25 and 50 p.s.i.a. runs was determined by the 
occurrence of critical flow. The heat transfer coefficients 


calculated from the data are estimated to have an accuracy of 


5%-10%, except for runs at very high flow rates w here the 
film temperature drop (7— 7) was low. For these latter data 
an accuracy of 15%-20% is estimated, excluding possibly two 
runs where the film temperature differences were as low as 
3.6°F, Generally the differences were greater than about 6.0°F. 

To determine if the data were in the convection regime, 
log-log plots of heat flux versus the film temperature difference 
(1 —T,) were prepared as in Figures 3 and 4. 
these lines is 1.0 at the higher flow rates, indicating that the 
heat transfer coefficient is independent of the heat flux. At 
lower flow rates the slope was found to increase with a decrease 
in either quality or flow rate, as in Figure 5. It is probable 
that this increase of coefficient with an increase in heat flux is 
due to a nucleate boiling contribution, for a convection mechanism 
does explain this phenomenon. At the higher flow rates nucleate 
boiling is probably suppressed because the wall temperature is 
too low to foster nucleation and bubble grow th, and it was found 
that Dengler’s correlation (8) for the minimum temperature 
difference required to initiate nucleate boiling, would predict 
this result. 

AT; = 10(G,V,/3600R,)°3 


At each pressure investigated, the data for the convection- 
controlled runs, smoothed by the above plots, were plotted on 
log-log plots of the heat transfer coefficient versus the superficial 

vapor mass velocity, as illustrated in Figure 6. The data were 
then found to fit equations of the form: 


herp = BCG, ooo c cc cccees. an 


where the exponent, 0.87 was the average slope for the four 
pressures studied, and the values of the pressure dependent 
constant, B, are listed in ‘Table 1. 


*A complete tabulation of data and results, as well as further details of 
the apparatus, procedure, and other phases of the study are given by 
Davis , 
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A maximum heat transfer coefficient was found to occur for 
each pressure and flow rate in the region of 80%-90% steam, 
but the phenomenon was not studied in sufficient detail to 
permit prediction of the maximum coefficient, because of the 


danger of burnout of the test section. Visual observation of 


the flow patterns indicated that for runs below the maximum 
heat transfer coefficient a thin liquid layer existed on the wall, 
and either annular flow or dispersed- annular flow occurred. 
The flow patterns were in general agreement with Baker’s 
flow patterns correlation‘?*), as shown in Figure 7. 


Correlation of Results 


Based upon the observed flow patterns and the correlation 
of the experimental data by the superficial vapor mass velocity, 
two models are postulated to predict the two-phase heat transfer 
coefficients for a vapor-liquid system. 


Separated-Annular Flow Model: I wo alternate approaches 

were developed based upon a separated-annular flow model - 

a slip ratio equation and a void fraction equation. For both 
jemand it 1s assumed that: 


(1) Nucleate boiling does not occur, because of suppression by 
two-phase forced convection heat transfer. 


(2) The resistance to heat transfer occurs within the thin layer 
of liquid at the wall. 


(3) The Reynolds number of the liquid and the physical prop- 
erties of the liquid determine the convection heat transfer 
rate from the wall to the liquid film. 

(4) An equation of the following form can be used to predict 
the heat transfer coefficients for heat transfer from the wall 
to the liquid film: 


hD, D.Gi\’ a) 
L 


.(2) 
ky ML k 


The liquid mass velocity in Equation (2) can be related to 
the vapor mass velocity by the equation of continuity for each 
phase and by the slip ratio definition: 


Pu, = G,, and pyu, = G_, 


a = u,/uz = slip ratio 
Solving for G;: 
G,px 
Got 
ap, 


Substituting in Equation (2) for G,: 


hD, DGa.\ (Cu 
= of PL ob .(2a) 
ki Miap, k L 


Because in the high vapor mass fraction range the vapor phase 
occupies practically all of the cross sectional area of the duct, 
the actual vapor mass velocity can be approximated by the 
superficial vapor mass velocity, and Equation (2a) can be written: 


AD, D.G,'piN’ "(& 


=a 
ki Miap, k 


.(3a) 
Because slip ratio data are not available for the range of pressures, 
qualities, and flow rates studied in the present investigation, it is 
not possible to assess the value of Equation (3a) per se for 
correlating the experimental data. Slip ratio data can be obtained 
from void fraction data to obtain the void fraction equation, as 
discussed below, but an alternate approach based upon Equation 
(3a) was developed to provide an empirical correlation of the 
experimental data. 

By including the values of the physical properties for each 
pressure studied in Equation (3a) and equating Equations (3a) 
and (1), the numerical values of a/a°-*? can be calculated for 
each pressure, as tabulated in Table 2. In the computations the 
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TABLE 2 


THE VARIATION OF a/a®:87 WITH HP RESSURE 


a/a®-87x108 


Pressure, p.s.i.a. 


1.051 25 
1.530 50 
2.385 100 
2.745 150 


equivalent diameter, D,, was taken as defined for single-phase 
flow: 


D 4c ross sec tional area for flow ) 
= ase 
"wetted perimeter 


However, there is no theoretical justification for this usage. 

The term, a/a°-’?, was found to be correlated by a log- log 
plot of the term versus the absolute pressure and also versus the 
liquid-vapor density ratio, p,/p,. The latter log-log plot results 
in an excellent fit of the variables by a straight line having a 
slope of 0.59, and this correlation can be included in Equation 
(3a) to give: 


hD, 0.59 DG , 0.87 C 0.4 
; = 0.060 ( % See ee c.... ee 


L Ps MiP, k L ; 


By relating the superficial vapor mass velocity to the total mass 
velocity, the equation can be written in final form: 
AD, pL 0.28 D.xG, 0.87 Cou o4 


= (0.060 - ong (4a) 
ky Ps Mr k L 


Equation (4a) was used to correlate the data in the convec- 
tion-controlled heat transfer region of the present investigation 
to within +15%, as shown in Figure 8. 

Although the range of values for the physical properties of 
the liquid ‘film was too small in the present study to permit 
independent determination of the effects of these properties on 
the heat transfer coefficient, inclusion of these properties in 
Equation (4a) produced a slightly better correlation of results 
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than was obtained from a generalized correlation of either a 
form similar to Equation (1), or to Equation (4a) used without 
the Prandtl number. 

The two-phase convection heat transfer data of Kvamme 
and Yamazaki“ and Dengler“:® for circular tubes were also 
correlated adequately by Equation (4a), as shown in Figure 9. 
(Because of the large scatter in the indiv idual point heat transfer 
coefficients measured along Kvamme’s test section, the average 
heat transfer coefficient for each of his runs was used.) Table 3 
lists the range of variables involved in the data correlated by 
Equation (4a). 

The ability of the correlating equation to predict the various 
data adequately indicates that the data are in agreement, that 
the equivalent diameter used is satisfactory, and that the model 
might be correct. However, because slip ratio data are not 
available, the model cannot be thoroughly tested, and the corre- 
lation must be considered as a purely ‘empirical approach, 

The separated-annular flow model can also be developed 
using void fraction data to calculate the slip ratio. The liquid 
void fraction is defined as: 

Rr = Azt/At = Axt/(AL + A, 


By use of equations of continuity for each phase, and the 
definitions of the void fraction and slip ratio, the slip ratio 
can be obtained in terms of the void fraction: 
x PL Ri 
(1—x) p, (1—Rz) 
Substituting in Equation (2) for a, the following correlating 
equation 1s obtained: 
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Figure 9—Correlation of the experimental data of Dengler, 


Kvamme, and Yamazaki. 
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Heat Flux Mass Velocity Equiv. Dia. 
B.t.u./hr.-ft.? lb. /hr.-ft.? ft. 
60,000-260,000 | 75,000-600,000 0.0324 
25,600- 50,200 85,000-189,000 0.0313 
27,600-198,000 41,000-284,500 0.0833 


In Equation (5), the void fraction equation, the exponents and 
the equivalent diameter are taken as in Equation (3a). The 
constant, 4, was determined to be 0.017 be fitting the equation 
to the data of the present study, using the liquid void fraction 
data of Smith and Hoe‘®, who correlated their data by the 
Lockhart-Martinelli method. Figure 10 is a comparison 
between the experimental data at a pressure of 50 p.s.i.a. and 
Equation (5). The correlation is unsatisfactory at both high 
and low qualities, but it cannot be determined w hether the heat 
transfer equation, Equation (5), is unsuitable or whether the 
void fraction correlation of Smith and Hoe is unsatisfactory in 
the high flow rate region. 

Equation (5) is similar to the correlating equation proposed 
by Kvamme®; 

hrp/h, = 1.5 [(1—x)/Rz]°% 


k DG, 0.F é; o.4 
ht = 0.023 — oe 


However, Equation (5) predicts two-phase heat transfer coefhi- 
cients that are about 20% higher than Kvamme’s correlation. 


where: 


Homogeneous Model: [he second correlating approach, 
the homogeneous model, was used because the flow pattern was 
aren or dispersed-annular for many of the heat transfer 
runs. With this model it is assumed that the heat transfer 
occurs to a very thin film of liquid on the wall, and the flow 
characteristics of the film are a function of the Reynolds number 
of the homogenous core. The equation used is: 


hD, D.pretre — Cou ai 
, eacaile 
ky Mre k L 


By considering a unit volume of duct, assuming no slip occurs, 
and that thermodynamic equilibrium exists between the two 
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Figure 10—Comparison of the void fraction model (Equa- 
tion 5) with the experimental data at 50 p.s.i.a. 
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Figure 11—Correlation of the experimental data at 50 
p.s.i.a. by the homogeneous model. 


phases, the following relations can be derived for the homo- 
geneous average linear velocity and the density of the homo- 
geneous mixture: 


urp = [x Pp + (1—x)/pi]G.... ein ical 
1/prp = x/p, + (1—x)/pr..... ee ..(8) 


It is assumed that the viscosity of the two-phase dispersed core 
can be defined in a manner analogous to the density. McAdams 
and Woods“® have used a similar definition in their two-phase 
pressure drop correlation: 


1 Mre = X/pP, (1 ¥)/ Pe . (9) 
The correlating equation developed from this model is: 


hD, DB LC, 
. = 0033 ( m 
L Mirp k L 


6.4 


Figures 11 and 12 show Equation (10) compared to the experi- 
mental data of the present investigation; the data are correlated 
to within +15%. bach curve on the figures represents a range 
of qualitics from 30%, to 90%, at some constant total mass 
velocity. 

Equation (10) is strikingly suniles to a correlation recently 
published by Groothuis and Hendal®® for heat transfer to air- 
water mixtures flowing in a circular duct 


Nu 0.029 Rep®*? Pr! ay / py)?! ; (11) 
Dp,u,' Dpuz' 
4 


where: Re: 
Me Hi 


Equation (11) can be rewritten as follows: 


hD. [= DG,\%" ( e (“ 0.16 
0.029 + (12) 
ky BK, i k Ma 


but: G,’ = xG,, and G,' = (1—x)G, 


Therefore: 


AD, D6." Cas" u\ 
= 9.029 ; (13) 
ky Mis k Me 


Iq uations (1Q) and (13) predict heat transfer coefficients within 
10%, of each other for the steam-water data taken in the present 
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Figure 12—Correlation of the experimental data at 100 
p.s.i.a. by the homogeneous model. 


investigation. ‘The agreement between the two equations is 
excellent, and indicates remarkable agreement between steam- 
water and air-water data for the wide range of flow rates, vapor 
mass fractions, and pressures studied. However, the highly 
empirical form of the equation does not attest to the theoretical 
correctness of the homogeneous model. 


Pressure Drop Correlation 


The pressure drop data taken with the rectangular duct 
of the present investigation were adequately correlated by 
the Lockhart-Martinelli correlation for isothermal two- phase 
flow®®, as shown in Figure 13. In preparing the experimental 
data for Figure 13, the single- phase friction factors needed were 
taken from a plot prepared from single-phase pressure drop 
measurements on water and superheated steam in the test section 
of this study; and the pressure drop due to acceleration was 
neglected. However, the acceleration pressure drop was esti- 
mated to be significant at higher heat fluxes and lower flow 
rates, where the quality increase in the test section was greatest. 
For most of the runs, because of the short heated length (6.0-in.) 
and short pressure drop length (4.0 to 11.0-in.), the ‘acceleration 
pressure drop was usually less than 15% of the total pressure 
drop (frictional plus acccleration). ‘The acceleration pressure 
drop was estimated by a momentum balance around the pressure 
drop section as suggested by Martinelli and Nelson®”, using 
both a separated- annular flow model and a homogeneous flow 
model. 


Conclusions 


‘The heat transfer coefficients for two-phase vapor- liquid 
flow at high vapor qualities can be predicted from equations 
similar to the Dittus-Boelter or Sieder- ate equations for single- 
phase flow. ‘Two models, a separated-annular flow model and 
a homogeneous model, have been used to develop such equations 
for two-phase flow. A slip ratio approach provided a correlation 
for much of the available steam-water data in the region of high 
vapor mass velocities: : 


AD, pr 0.28 dD, xG; 0.87 ran 04 
0.060 a (4a) 
ky Pe Mi k L 


Iquation (4a) has been used to correlate the data of the pres- 
cnt investigation to within +159, and the data of Kvamme, 
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Figure 13—Experimental two-phase pressure drop data. 


Yamazaki, and Dengler to +20%. An alternate equation based 
on a void fraction “approach (Equation 15) proved relatively 
unsatisfactory. A homogeneous model correlated the data to 
approximately +15%, and the equation developed and given 
below was found to be in excellent agreement with the Groothuis- 
Hendal correlation for air-water flow: 


\ 0.87 


hD, D.G é 
— = 0.033 : a .. (10) 
ki Mrp 


The non-isothermal two-phase pressure drop data were found 
to be correlated adequately by the Lockhart-Martinelli correla- 
tion for isothermal flow. 
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Nomenclature 

a = [Dimensionless constant 

A = Cross sectional area of tube or duct, ft.? 

b = Dimensionless exponent 

3 = Pressure dependent constant in the equation, hyp = 
B(G,')°:87 

¢ = Dimensionless exponent 

C, = Heat capacity at constant pressure, B.t.u./Ib.-°F. 

D = ‘Tube diameter, ft. 

D, = Equivalent diameter of a tube of duct, ft. 

G Mass velocity, Ib./hr.-ft.? 

G’ Superficial mass velocity, Ib./hr.-ft.? 

h Heat transfer coefficient, B.t.u./hr.-ft.2-°F. 

hy Heat transfer coefficient based on the total flow rate as 
liquid, B.t.u./hr.-ft.2-°F. 

hrp ‘Two-phase heat transfer coefficient, B.t.u./hr.-ft.2-°F. 

k Thermal Conductivity, B.t.u./hr.-ft.2-°F. /ft. 

L Superficial mass velocity of liquid phase, Ib./hr.-ft.2 (In 


Baker correlation, Figure 7) 


R, A,/A, Liquid void fraction, dimensionless 
R, A,/A, Vapor void fraction, dimensionless 
Res Reynolds number defined by Groothuis and Hendal, 
dimensionless 
l Pemperature, °F. 
AT; The minimum temperature difference required to initiate 
nucleate boiling, °F. 
u Average linear velocity of a phase, ft./see. 
upp Average linear veloc ity ola homogeneous mixture, {t./sec. 
V, Specific volume of the liquid phase, {t.4/Ib. 
Mass flow rate, Ib./see. or Ib. /hr. 
\ Vapor mass fraction, dimensionless 
(AP/AL) rp ; 
X, Martinelli’s pressure drop parameter, 
(AP/AL), dimensionless 
Greek Letters 
a The slip ratio, dimensionless, equals u,/uy 
\ Baker’s flow pattern parameter, dime nsionless 
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bu = Viscosity, lb./hr.-ft. 


Mrp = Viscosity of a two-phase homogeneous mixture, Ib. /hr.-ft. 
p = Density, lb./ft.3 
Pre = Density of a two-phase homogeneous mixture, lb./ft.3 


(AP, AL)rp oe 


gd. = aie Martinelli’s pressure drop parameter, 
(AP/AL)z dimensionless 

Y = Dimensionless flow pattern parameter used in Baker 
correlation 

Subscripts 

b = Refers to a physical property at the bulk temperature 

g = Refers to the vapor phase 

L = Refers to the liquid phase 

s = Refers to saturation temperature 

t = Refers to a total value 

TP = Refers to two-phase flow or to a homogeneous two-phase 
property 

w = Refers toa physical property at the wall temperature or 
to the wall temperature 
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Fluid Friction and Heat Transfer in 
Cylindrical Pipes: Relationship between 
Lumped and Distributed Parameters’ 


JULIAN C. SMITH? 


Physical phenomena may be described from the 
“microscopic” point of view, in terms of distributed 
parameters, or from the “macroscopic” viewpoint, 
using lumped parameters. Engineering correlations 
almost always involve lumped parameters. A lumped 
parameter is equal or proportional to some average 
value of the corresponding distributed parameter. 
Like any average quantity it gives no information 
regarding the form of the distribution. The use of 
lumped parameters therefore requires some assump- 
tions, expressed or implied, regarding the form of 
the distribution within the system. 

The characteristic lumped parameters used in 
problems of heat transfer and fluid friction are over- 
all or average fluxes of heat energy or momentum, 
arbitrarily defined as described in this paper. Com- 
mon dimensionless groups, such as f, Ny., Ns, and 
others are ratios of these energy or momentum fluxes. 
Such dimensionless groups, since they are ratios of 
lumped parameters, are themselves lumped para- 
meters. In their definition, therefore, several assump- 
tions regarding the distributions of velocity and 
temperature are implicit. When these assumptions 
are not valid, correction factors such as diffusivity 
ratios or length-to-length ratios are needed to allow 
for deviations from the distributions used for refer- 
ence. 

Using this approach the form of empirical equa- 
tions becomes easy to predict. Dimensional analysis 
is not needed. In addition, this approach suggests 
that some accepted correlations may contain weak- 
nesses not predicted by other methods. 


ene occurring in a physical system may be analyzed 
in terms of the behavior of the entire system, or in terms 
of events taking place at individual locations within the systems. 
From the internal or “microscopic” viewpoint transfer phenom 
ena are described by differential equations relating the distributed 
parameters of the system. From the external or “macroscopic” 
point of view phenomena are described by lumped parameters 
which are characteristic of the whole system. Lumped parame- 
ters are usually restricted in their application to a particular 
system or to a small group of highly similar systems. However, 


iManuscript received September 24, 1960; accepted March 6, 1961 
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ing, Cornell University, Ithaca, New York, U.S.A 
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they are of direct practical use, whereas distributed parameters 
ordinarily are not. Engineers are therefore much more often 
concerned with lumped parameters than with distributed para- 
meters. 

Full understanding of external behavior, however, demands 
a knowledge of internal events and of the relationship of one to 
the other. The transition from microscopic to macroscopic, 
from distributed parameters to lumped parameters, may be 
made in several ways. In a few cases the differential equations 
may be integrated analytically. Usually, however, this is not 
possible and the transition must be made by dimensional analysis 
of the pertinent differential equations” or by application of the 
principles of similitude®:®, Unfortunately these methods do 
not make it easy to visualize the underlying physical relation- 
ships, nor do they emphasize the assumptions involved. This 
paper suggests a way of interpreting some commonly used 
lumped parameters in terms of the corresponding distributed 
parameters, bring out these relationships and assumptions. 


Fluid Friction 


‘The relationships between lumped and distributed parameters 
in a flowing fluid will be developed for a specific and rather 
simple case that of, steady fully developed flow of an in- 
compressible fluid in a horizontal cylindrical pipe of constant 
cross section. In this system the relationships are reasonably 
clear. Many of the commonly used lumped parameters appear 
to have been originally defined with reference to this system, 
though they have since been applied with more or less 
success to other more complex situations. ‘The concepts 
presented here may also be applied, with some modifications, 
to other physical systems. 


Distributed Parameters. Consider the cross section of a 
horizontal cylindrical pipe of diameter D in which an incom- 
pressible fluid of density p is flowing in steady, fully developed 
flow. Assume for the present that viscosity ws is constant. 
Flow may be laminar or turbulent, but the bulk motion of the 
fluid is unidirectional and parallel to the pipe axis. 


The distributed parameters are the fluid velocity u, pressure 


p, and shear stress 7, which is related to the transfer flux* of 


momentum ju by the equation 
ju = Reka. (1) 


Since # and p are assumed constant they need not be con- 
sidered as distributed parameters. Furthermore, p is nearly 
constant across the pipe cross section. Only the velocity and 
the transfer flux of momentum vary significantly from point to 
point. 


°A “transfer flux’’ is the flux resulting from a potential gradient, A 
“convective flux’’ is the flux resulting from motion of the fluid. 
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Figure 1 —Profiles of velocity and momentum flux in 
turbulent flow. 


The nature of this variation is indicated by the velocity and 
flux profiles in Figure 1. The velocity is zero at the pipe wall 
and a maximum at centerline. Flux jx is zero at the centerline 
and rises to a maximum of j;» at the pipe wall. If p is constant 
across the cross section the variation of jy is linear with radial 
distance from the pipe axis. 


In general the total momentum flux is made up of two parts: 
the flux caused by viscous action, jyy, and the flux caused by 
turbulent action, juz. When flow is laminar jwr is zero and 
juv = ju. With turbulent flow the flux by viscous action is 
normally a very small part of the total flux except near the pipe 
wall. Nevertheless it is not zero except at the pipe axis. It 
ty pically varies with radial distance as indicated in Figure Fr. 
The viscous flux juv is therefore a distributed parameter, given 
by the equation 


Sadie Se re os Nis ee een her (2) 


Quantity jue, a particular local value of a distributed par- 
ameter, 1S of especial interest to the engineer, because from it 
he can find the frictional pressure drop by the relation 


1 4Ty tiuw 
i a (3) 
dx D Deg, 


where 7, is the shear stress at the pipe wall. When the velocity 
profile may be calculated mathematically, as in the fully devel- 
oped isothermal laminar flow of a Newtonian fluid in a cy lindrical 
pipe, T, may be calculated directly from the relation 


du 
Twhe = w bine soennacwee) 
dy 
oy. 
where (du/dy)» is the velocity gradient at the wall. This leads 
to the Hagen-Poiseuille equation for frictional pressure drop in 
laminar flow. When the shape of the velocity profile is unknown, 
however, (du/dy)» cannot be found directly, and empirical 
correlations using lumped parameters must be employed. 


Lumped Parameters in Fluid Flow. Ay important: group 
of lumped parameters is comprised of the overall convective 
Huxes; that is, the convective flux caused by the motion of all 
the fluid passing through a pipe cross section. First is the average 
velocity P, which is the overall convective flux of volume (total 
volume flowing through a unit cross sectional area per unit time). 
Other overall convective fluxes are the products of P with 

“concentration”. The overall convective flux of mass is /p or 
G, since p is the “concentration” of mass. The average ‘con 
centration” of momentum in an incompressible fluid equals ol, 
and the overall convective flux of momentum, J, is given by* 


Jy = pvr : pV? (5) 


. o yy 2 “ 3 
*Rigorously Ja = pu?* or pV? B, where PB is a “convective flux 
correction factor.”’” For fully developed laminar flow in a cylin 
drical pipe 8 = 0.75; for turbulent flow it is almost unity. 
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Figure 2—Definition of “average” velocity gradient. 


An important but loosely defined lumped parameter is the 
average value of j, in the pipe. As shown later, this average 
is not used without modification and thus for present purposes 
it may be defined in any convenient way so long as the same 
definition is alw ays used. One such way ts illustrated by Figure 
2, which shows a typical velocity profile i in the fluid. Point P 
is the point on the ‘profile where the local velocity u equals v. 
The slope of the line OP may be taken as a measure of the 
average velocity gradient in the fluid, and an “‘av erage”’ value 
of jury defined by 

- V , 
Juv = io Ay . (6) 
where Ay is the radial distance from Point P to the pipe wall. 

This “‘average’’, however, is valueless as a lumped parameter, 
for Ay is oulinaseits unless the shape of the velocity profile is 
known. A useful lumped parameter, however, is formed by 


defining a ‘‘modified average”’ transfer flux, 7 uy, by the equation 


* =, my I 


-= = — 7) 
Jy Jmy D BOD (/ 


where D is the pipe diameter. Flux j vy May be regarded as 
the viscous flux of momentum corresponding to a “modified 
average’’ velocity gradient, which is the slope of line OQ in 
Figure 2. 

Flux Ratios. Friction Factor. \ secondary or derived 
lumped parameter is the ratio of the scalar magnitude of the 
transfer flux at the wall to that of the total convective flux of 
momentum. These flux vectors are shown in Figure 3. This 
ratio, by definition, ts proportional to a friction factor. Specifi- 
cally, the F anning friction factor is defined by the equation 


Mw 
» J 
é /M 
Substitution from Equations (3) and (5) into Equation (8) 
gives the more familiar definition 


dp Deg 
m 7 (9) 
dx 2p? 
| 
-* Ju 
ong 
Yiny 
| PIPE WALL 
Ju, 


Figure 3—Momentum flux vectors in pipe flow. 
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Ratio of Convective Flux to “Modified Average”’ Transfer 
Flux. Pipe Reynolds Number. Another secondary lumped 
parameter is the pipe Reynolds number, which is proportional 
to the ratio of the total convective flux to the “modified average’ 
transfer flux. It is defined by 


Substitution from Equations (5) and (7) into Equation (10) 
gives 


Pe _ DP, 
a Be. ita ae (11) 
pV /D Ke 


Note that Juv contains an average flux jas which, however 
defined, may be considered proportional to the true average 
momentum flux in the fluid, and also a term Ay/D which depends 
on the shape of the velocity profile. Thus the significance of the 
Reynolds number depends on the shape of the velocity profile, 
even when all other characteristics of the system are held 
constant. As stated by McAdams: “In two geometrically 
and kinematically similar systems, equal Reynolds numbers assure 
dynamic similarity.” This assumption of kinematic similarity 
in the use of the Rey nolds number is often overlooked. 


Relationship between Friction Factor and Reynolds 
Number. The friction factor is proportional to the lumped 
parameter J,, and the local value of the distributed parameter, 
the transfer flux at the wall, jv. The Reynolds number is also 
proportional to Jy and to another lumped parameter, 7 ‘sy, 
which is a ‘“‘modified average” value of the viscous transfer 
flux of momentum. If it is assumed that the transfer flux at 
the wall depends on the “average” viscous transfer flux, then 
there must be a functional relationship between f and Nx. 

This relationship is simple in fully developed laminar flow, 
in which the transfer flux everywhere in the fluid is by viscous 
action and the velocity profile is parabolic. The conditions of 
kinematic similarity are met, and there can be no variation in 
the ratio Ay/D. Under these conditions the Reynolds number 
and friction factor are both proportional (though in different 
ways) to the same physical quantities. The one-to-one relation- 
ship between them is given by the simple equation 

f= ad (12) 
Nace 2 

When flow is turbulent the relationship is more tenuous, 
for f and Nx, are no longer measures of the same physical 
quantities. Unless the average viscous flux always bears the 
same relation to the total transfer flux, the Reynolds number 
cannot be uniquely related to the friction factor. This require- 
ment is satisfied if kinematic similarity is assumed. Anything 
which changes the shape of the velocity profile, however, 
changes the friction factor even though the Reynolds number 
is held constant. 

The roughness of the pipe wall has this effect and the con- 
ventional friction-factor plot contains lines for different pipe 
roughnesses. For a pipe of given roughness a unique relationship 
between f and Nx, is assumed; that is, it is assumed that at a 
given Reynolds number the shape of the velocity profile 1 is always 
the same. This may or may not be a valid assumption. It is 
equivalent to assuming that ‘the Prandtl mixing length is inde- 
pendent of molecular viscosity, or, alternatively, that the eddy 
diffusivity of momentum, at a given point in the fluid and at a 
given Reynolds number, bears a constant relationship to the 
molecular viscosity. If, for example, the fluid viscosity is in- 
creased but Nx. is held constant, by changing D or V, then ey 
at corresponding points in the fluid must increase in exact 
proportion to the increase in v. Otherwise the shape of the 
velocity profile would change and f would not bear the same 
relationship to Ny. 

Exactly proportional increase of the turbulent and molecular 
diffusivities at all points seems improbable, and some change 
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in the shape of the velocity profile with molecular viscosity 
might be expected. The possibility of such a change is suggested 
in the analysis of turbulent-flow velocity profiles by Hanratty 
and Flint“, 


Published data do not indicate any variation of f with v at 
a given N,.. The range of viscosities studied, however, has 
been rather narrow. From a practical standpoint the variation, 
if it exists, is probably insignificant. Furthermore, it is difficult 
to study, for to create turbulence in a highly viscous liquid 
causes so much energy loss by viscous dissipation that the 
assumptions of isothermal flow and constant viscosity are no 
longer valid. On the other hand, as pointed out by White and 
Churchill“, the fact that experimental measurements (especially 
of transfer rates) appear to agree with an integral equation is 
no assurance that the differential model which leads to the 
integral equation is valid. It seems probable that in the inter- 
mediate range of turbulence f and Ng are not uniquely related 
for a pipe of given wall roughness, but existing data are not 
sufficient to indicate whether or not this is true. 


Though unimportant in most practical engineering calcula- 
tions, the question is conceptually important. Many empirical 
correlations, in heat transfer as well as fluid dynamics, may be 
regarded as using the fully developed velocity profile of a 
constant-viscosity fluid as a reference or standard. A given 
Reynolds number under these conditions is assumed to imply a 
velocity profile of a = shape. Actually there are two 
“standard” profiles: one for laminar and one for turbulent flow. 
Any variations in the shape of the turbulent-flow profile with 
molecular viscosity are ignored. 

At high turbulence the viscous transfer flux is a negligible 
part of the total, and variations in its average value have no 
effect on the transfer flux at the pipe wall. Changes in j ‘yy 
do not affect jarw. As shown by Equations (8) and (10), f and 
Nee in this region are therefore measures of unrelated physical 
quantities and are independent of each other. 


Effect of Temperature Variation. \When the temperature 
of a Newtonian fluid varies across the pipe, the changes in 
viscosity alter the form of the velocity profile. This effect is 
illustrated in Figure 4. The velocity gradient at the pipe wall 
becomes steeper, and the momentum transfer flux at the wall 
becomes 

du 


jmMw = — pe ee ae ‘ a it FS) 
dy i 


This change is reflected in the friction factor. The Reynolds 
number, however, does not change, except that the “modified 
average” transfer flux now involves still another lumped para- 
meter, the bulk average viscosity ff, defined as the viscosity 
at the bulk average temperature, 7. Actually the viscosity 1s 
now a distributed parameter, varying in a way dictated by the 
variation in fluid temperature. The Reynolds number, therefore, 
is less of a true measure of conditions in the fluid than with 
isothermal flow, for it makes no allowance for changes in either 
the velocity profile or the temperature profile. Furthermore, 
: T changes appreciably with distance along the pipe, # must 

“Jumped” along the pipe as well as across it. Ng. does not 
tt into account the form of the temperature change along the 
pipe: it remains a single “‘average”’ value. 

Since f changes while Nz, does not, a correction factor must 
be included in the equation relating them, to restore kinematic 
similarity. This Sieder-Tate correction factor is 1.1 (/pw)-?® 
for laminar flow and (f/u,)°- for turbulent flow™. It is 
assumed, in applying this factor, that the change in velocity pro- 
file is always the same function of this ratio of a lumped parameter 
fi to the local value of the distributed parameter wy. This is 
unlikely, especially for turbulent flow, and probably the correc- 
tion factor should also depend on the absolute magnitude of the 
fluid viscosity. Again, however, the effect would be difficult to 
establish experimentally. 


The Canadian Journal of Chemical Engineering, June, 1961 


2.25 


2.00 


1.75 


1.50 
u/V 


1.25 


1.00 


0.75 


0.5C 


0.2s 


Figure 


Iso 
of an it 
varies 
gradier 
from t 
viscous 


fluids, 


where 
imw, be 
flux is 


where 
overal 
numbe 


TI 
Metzr 
Nr. fe 


The 





2.25 


2.00 


1.75 


1.50 
u/V 


1.25 


1.00 


0.75 


0.50 


0.25 


0 0.2 0.4 0.6 0.8 1.0 


y/ tw 


Figure 4— Effect of temperature gradients on velocity 
profile. 


Isothermal Flow of a Non-Newtonian Fluid. [pn the flow 
of an incompressible non-Newtonian fluid the apparent viscosity 
varies across the pipe, since it is a function of the velocity 
gradient. The shape of the velocity profile therefore differs 
from that in a Newtonian fluid, as shown in Figure 5. The 
viscous transfer flux of momentum, for many non-Newtonian 
fluids, is given by an empirical equation of the form 


du 


juv Oh Ee 5 oe ath eo ic Sc chnals Oaks (14) 
dy 


where K and 7 are constants. The momentum flux at the wall, 
iuw, becomes — K(du/dy)", and the “modified average” viscous 
flux is 


» K V\" (dyv\" nih ‘i 
— iene aan =— tot oer 
J MV Ay D D tees ( 


where D/Ay i is the average velocity gradient in the fluid. The 
overall convective flux, Ji, remains ‘equal to pV’*. a he Reynolds 
number, then, by substitution i in Equation (10), i 

Ju D*P2 "p 


Nre = = — 
jmv K 


This is the Reynolds number Nx° used by Dodge and 
Metzner“, When flow is laminar the relation between f and 
Nee for a non-Newtonian fluid is 
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Figure 5— Velocity profiles in laminar flow of non- 
Newtonian fluids. 
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Function $(7) is analogous to the Sieder-Tate correction 
factor for the non-isothermal flow of a Newtonian fluid. It 
restores kinematic similarity by “correcting” the shape of the 
velocity profile to that obtained in the fully developed laminar 
flow of a Newtonian fluid. 

With — flow the shape of the velocity profile also 
depends on 7, although presumably the ‘“‘correction factor” 
would not ‘ the same as that for laminar flow. An equation 
relating f, Ne, and 7 has been proposed by Dodge and Metzner®. 
They assumed that the shape of the profile is ‘independent of K 
in the turbulent region; as with a Newtonian fluid (for which 
K = p) it seems probable that this is not quite true. 


Heat Transfer 


In heat transfer by conduction thermal energy is transferred 
by virtue of temperature gradient. Since temperature is a 
scalar the heat fluxes are simple vectors, not tensors as are 
momentum fluxes (stresses) in fluid dynamics. The situation 
in heat transfer thus appears simpler, at first glance, than in 
fluid flow, and the lumped parameters might appear to involve 
fewer assumptions. Actually the reverse is true. In heat transfer 
to or from a fluid flowing in a pipe, the transfer flux at wall is 
nearly always more dependent on the fluid motion, especially 
in turbulent flow, than on thermal conduction. Heat transfer 


109 








is thus influenced by two superimposed and interrelated fields: 
temperature gradient vectors and the stress tensors. 


Distributed Parameters. The important distributed par- 
ameters are the temperature and specific enthalpy of the fluid, 
the thermal diffusivity, and the transfer flux of thermal energy. 

Consider the fluid passing through a cross section of the 
pipe. Let the temperature at the pipe wall, 7,, be used as the 
datum for calculating the specific enthalpy H. Assume that 
there is no phase change in the fluid and that its heat capacity ¢c, 
is constant. Then 

r 
H= | edT =0(T — Ty)... (18) 


« 


The enthalpy per unit volume of fluid — the ‘‘concentration”’ 
of enthalpy — is pH or pc,(T — T.). The total transfer flux 
of thermal energy, jz, is related to the gradient of this con- 
centration when p is constant by 


d(pH) 


jn = — (a + €n) — 
i dy 


dT 
AE SO (19) 
7 


This flux varies from zero at the pipe axis to a maximum 
value, jye, at the pipe wall. If all the transfer is by conduction, 
€, is zero and at all points 

: dT b dT (20 
Hv = — Qply a . ) 
, dy dy 

The heat flux at the wall is related to the individual heat 

transfer coefficient / by the equation 


jaw = MT — Te)... (21) 


Lumped Parameters. An important lumped parameter is 
the average or “mixing-cup” temperature 7, which is related 
to the average enthalpy H of all the fluid in the cross section 
by the equation 

H =c,T — T.) (22) 


The average concentration of enthalpy is p/7, and the overall 
convective flux of enthalpy, J, is* 


Ju = Vo = Voc,(T —- T.) (23) 


_ The “modified average” transfer flux of thermal energy, 
J'nv, 1s defined in the same way as the corresponding momentum 
flux. It is 


ae k AT Ay 

Say * Ay . D 
=-k r= (24 
- D } 


where (T — T,) Ay is the “‘average’’ temperature gradient in 
the fluid indicated in Figure 6. 


Flux Ratios. Dimensionless Groups. [he ratios of the 
scalar magnitudes of the several heat fluxes, as thus defined, 
are proportional to common dimensionless groups. Thus from 
Fquations (21) and (23) 


jHw WT — Ty) h Vy (25 
* > T * = . - 9) 
Ju Vp Ai a €,G 
From Equations (23) and (24) 
*As with the convective flux of momentum, this is not quite 


true. The overall convective flux of enthalpy is rigorously equal 
to pc, V(T Tw) OF Ply (7 7.)/B, where B is the “con- 


vective flux correction factor”’ 
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Figure 6—Definition of “average”? temperature gradient. 


J Voc(T — Te)  DGe 
= aoe = peo f ae = ee = Np. a Grow We tales (26) 
juv k(IT — T,)/D k 





Also, 
Jie h(T — Ty) hD S 
oe et cen te ee Ni i 5c RD 
Juv k(T —Ty)/D k 

Relationships among Dimensionless Groups. The lumped 
parameters Ns, and Np, are related in the same way as the cor- 
responding parameters for fluid friction, f and Ne When the 
thermal transfer flux is essentially all by molecular action, jy» 


and Sor are uniquely related and a unique relationship should 
exist between Ns, and Np,, or, since Ny, = Ns; Np, between 
Ny, and Np,. This assumes that the shape of the temperature 
profile is “fully developed” and is always the same at a fixed 
value of Np. It is apparently a valid assumption for heat 
transfer to liquid metals, as shown by the empirical equation 
of Lyon 

Nyy = 7 + 0.025 Np,®'8.... sie o Sa 


When flow is laminar the heat transfer is essentially all by 
conduction, but the temperature profile continuously changes 
with length and never becomes fully developed. The relationship 
is then of the form 


; E 14D 
Nwu = (Nez) = O( Np,) aL a (29) 


where Ng, is the Graetz number, which equals (rD/4L) Np,. 
The quantity #D/4L may be regarded as a correction factor, 
lumped over the pipe length, which allows for the difference 
between the changing temperature profile and the steady ‘fully 
developed” profile used in the definition of Np. 


With turbulent flow the temperature profile may be con- 
sidered fully developed, provided the temperature change along 
the pipe is small. The shape of the profile, however, is dictated 
not only by the thermal conductivity of the fluid but also by 
the fluid motion, and is therefore influenced by the molecular 
and eddy viscosities. Distance Ay, which may be considered 
constant for fully developed isothermal fluid flow, is no longer 
constant. For this reason Np, which contains this variable Ay 
implicit in its definition, is not as generally useful as Nx. 


If the temperature and velocity profiles, an chance, have 
the same shape, so that at any point 1, ‘V = T/T, the distance 
Ay in the definition of Nx and Np, + hele (7), (10), (24) 
and (26)) is the same for both dimensionless groups. The 
relationships between the overall convective flux and the transfer 
flux at the pipe wall are also the same. This is to say 


Jitw - Jie (30) 
Jun Ju 
‘This means, from Equations (8) and (25), that 
,=f/2.... (31) 


Equation (31) is a statement of the Reynolds analogy, and is 
approximately true for heat transfer to gases. 
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Figure 7—Relationship of N,, to temperature and velocity 
profiles. 


With most fluids, however, the profiles of temperature and 
velocity do not have the same shape. Correlations among the 
dimensionless groups may then be considered to use the fully 
developed constant-viscosity velocity profile as a standard to 
which the t — profiles are “corrected.” It is assumed 
that at a given N,, the velocity profile always has the same shape. 

The shape of the temperature profile is “corrected” to that 
of the velocity profile by using the ratio of the molecular diffu- 
sivities of momentum and thermal energy. That is to say, it is 
assumed that 


Ayu v 
— = Gata dip ce bat erhnn eed ih as 32) 
Ayu ® a ( ; 


Ratio v/a is the Prandtl number, N;,. The assumption stated by 
Equation (32) is somewhat dubious, since v and a refer only to 
molecular transfer and ignore any effect of turbulence. F ‘urther- 
more, v and a as normally used are themselves lumped parame- 
ters, the “bulk average” values for the fluid, and do not consider 
how the distributed values actually vary across the cross section. 

\s shown by experiment, however, the shape of the temper- 
ature profile varies significantly with Np,"®. This is illustrated 
in Figure 7. Possibly the effect of other variables is too small 
to establish experimentally. Assuming Equation (32) is valid, 
Equation (31) 1s modified to 


Ns: W(Np-) = f/2... ... (33) 


For turbulent flow, as found empirically, ~(Np,) = Np,*, 
and 


Ns Np,?!8 = f/ Res .(34) 


Since the temperature is not constant the shape of the velocity 
profile must also be ‘‘corrected’’ to allow for changes from the 


constant-viscosity profile by using the ratio Z/p,, and 


f 
Ns Np? = ° 
2 \ Be 
or 
bu 0.14 f 
Ns: Np?" — =jo=*.. (35) 
vi 2 


where jc is the Colburn “‘j factor’. Equation (35) is a statement 
of the Colburn analogy. Substitution of empirical relations 
between f and Nx, gives equations relating / and Nx. 
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The form of empirical correlations for turbulent-flow heat 
transfer may also be established as follows. As previously 
discussed, if the shape of the temperature profile depended only 
on Np,, then Ny, would be a unique function of Np.. Since 
this is not true, it is necessary to “‘correct’’ the shape of the 
temperature profile to that of the reference velocity profile. 
The correction factor is assumed to be a function of Np,. Then 


Nyw = X(N, Ne,)...-. (36) 
Since Np. = Nae Np,, Equation (36) becomes 
Nwu = X(Nre, Ney) 6h 


A recent empirical equation in this form is that of Friend 
and Metzner“ 


Nyu = 0.022 Ng,°-8 Np,?-#....... . (38) 


Summary and Conclusions 


(1) Heat transfer and fluid friction, even in fully dev eloped 
flow in a cylindrical pipe, involve complex interrelationships 
among a number of distributed parameters. Engineering correla- 
tions are almost always based on lumped parameters. Lumped 
and distributed parameters for simple cases of heat transfer and 
fluid friction in a cylindrical pipe are listed in Table 1 

(2) A lumped parameter is equal or proportional to some 
average value of the corresponding distributed parameter, and 
like any average gives no indication of the form of the distribu- 
tion. Ifa given lumped parameter depends on the form of the 
distribution, some assumption as to the form of the distribution 
must be made in defining the lumped parameter. 

(3) Common dimensionless groups f. Nae Nee Nsw 
etc. — are lumped parameters. They may be regarded as being 
proportional to the ratios of the scalar magnitudes of arbitrarily 
defined flux vectors. Np, for heat transfer is the analogue of 
Ne. for momentum transfer. 

(4) Nee implicitly depends on the form of the velocity distri- 
bution, and its use implies kinematic similarity. Np, and Ny, 
similarly depend on the form of the temperature distribution. 


(5) When kinematic similarity is assumed, consideration of 
the flux vectors leads directly to the form of empirical equations 
relating f and Nz,. When the shape of the temperature profile 
is assumed constant, the same reasoning leads to the form of 
equations relating Ns; and Np. 


TABLE 1 


LUMPED AND DISTRIBUTED PARAMETERS IN SIMPLE CASES 
oF FLurp FRICTION AND HEAT TRANSFER 


a Important 
Assumed | Distributed _— 
Local 


Constant | Parameters , 
Values 


Lumped 
Parameters 


\. Flow of constant-viscosity incompressible fluid 


v u Jin I = 
p JM G=Vp 
Jy = VG = pP? 
ren 
f Jie Jv 
| Nr Ju /F uy 


B. Heat transfer in forced convection of incompressible constant- 
viscosity fluid. 
| | 


v r jie r bis 
» A 
a in Ju = Volt = Voc lT—Te) 
Np, = v/a Juv jay 
Vs | tw din 
Vv; Ju jm 
| Vy, lHw/ J WN 


Ill 








(6) For situations in which the shape of the velocity profile 
or temperature profile cannot be considered constant, empirical 
correlations may be considered as being based on standard or 
reference profiles to which the actual profiles are ‘‘corrected”’ 


(7) The standard profiles are the velocity profiles found 
with steady fully developed flow of a constant-viscosity fluid. 
Two standards are used: one for laminar and one for turbulent 
flow. It is assumed that in turbulent flow the shape of the 
velocity profile depends only on Nx,, although this assumption 
appears somewhat dubious. 

(8) The “‘correction factor’’ which restores kinematic simi- 
larity in the non-isothermal flow of a Newtonian fluid is the 
Sieder-Tate correction factor. That for the isothermal laminar 
flow of a non-Newtonian fluid is @(7). 

(9) The Prandtl number performs a similar function in 


equations for turbulent heat transfer, “correcting” the shape of 


the temperature profile to that of the reference velocity profile. 


(10) Expressions for heat transfer in laminar flow involve 
New which may be regarded as the product of N,, with a 
“correction factor”, #D/4L. This factor allows for the differ- 
ence between the shape of the developing temperature profile 
and that of the reference profile. 

(11) Although adequate for most engineering purposes, some 
“correction factors” appear to be incomplete measures of changes 
in profile shape. Empirical correlations should perhaps contain 
other functions involving fluid properties, but such secondary 
“correction factors” would be difficult to establish experiment- 

ally. More desirable would be better analytical solutions to 
problems in turbulent transfer. ‘Traditional empirical correlations 
are based on oversimplifications of highly complex phenomena, 
and are therefore applicable only to relatively simple systems. 
This is not to deny their usefulness nor their importance, for 
they provide solutions to problems which can be solved in no 
other way. In the longer view, however, they are merely 
useful expedients which should be discarded when more powerful 
generalization, have been developed. They are but w ay-stations 
along the road to a better understanding of turbulent transfer 
processes. 


Nomenclature 
, = Heat capacity at constant pressure, B.t.u./Ib., °F. 

D = Pipe diameter, bE 

f = Fanning friction factor 

G = Mass velocity, lb./sec., ft.? 

gc = Newton's-law aiaialee factor, lb.-ft./lb.y-sec.? a 

H = Specific enthalpy above 7, as datum, B.t.u./lb.; (7, 
bulk average value 

h = Individual heat transfer coefficient, B.t.u./hr., ft.2, °F. 

J = Overall convective flux; J, flux of enthalpy, B.t.u./sec., 
{t.2; Jy, flux of momentum, Ib./ft., sec.? 

i = Transfer flux; jy, by molecular action; jy, flux of thermal 


energy, B.t.u./sec., {t.?; 7’”7y ‘modified average’ flux of 
thermal energy, defined by Equation (24); j4, momen- 
tum flux, Ib./ft., sec.2; jyy, ‘average’? momentum flux, 
defined by Equation (6); 7’yy, ‘‘modified average”’ 
momentum flux, defined by Equation (7). 
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ic = Colburn “‘j factor”’ 
K = Coefficient of non-Newtonian fluid, (Equation (14)) 
k = Thermal conductivity, B.t.u./sec., ft.2, (°F./ft.) 
L = Pipe length, ft. 
Ne: = Graetz number 
Ny, = Nusselt number 
pe = Peclet number 
pp = Prandtl number 
Nre = Reynolds number 
Ns = Stanton number 
n = Coefficient of non-Newtonian fluid, (Equation (14)) 
p = Pressure, lb.,/ft.2 _ 
r = Temperature, °F.; 7, bulk average temperature 
u = Local fluid velocity, ft./sec. 
x = Axial distance 
y = Radial distance 


Greek Letters 


a = Thermal diffusivity, ft.2/sec. 

8 = Convective flux correction factor 

Ay = Distance used in definition of Juv (Equation (6)) 

= Eddy diffusivity, ft.2/sec.; €”, of thermal energy; €,, of 
momentum 

bu = Absolute viscosity, lb./ft., sec.; @, bulk average vis- 
cosity 

v = Kinematic viscosity, ft.?/sec. 

p = Density, lb./ft.’ 

T = Shear stress, lb.,/ft.? 


¢,x,¥ = “Functions of” 


Subscripts 

M = For momentum transfer 

H = For heat transfer 

7 = By turbulent action 

V = By molecular or viscous action 
w = Value at pipe wall 
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Fundamental Aspects of Solids‘Gas Flow 


Part VI: Multiparticle Behavior in Turbulent Fluids’ 


L. B. TOROBIN? 


At Reynolds numbers approaching unity the par- 
ticles in a solids-gas system acquire a turbulent 
motion in response to the turbulence of the ambient 
fluid. This regime is of particular importance and a 
discussion of recent theoretical and experimental 
attempts to describe it is presented. The reaction of 
the fluid phase to the presence of the particles is 
equally important. The extent to which the mean and 
turbulent properties of the fluid motion are affected 
by the presence of the solids phase is governed by 
the particle concentration and the fluid and solids 
Reynolds numbers. The experimental evidence is as 
yet inconclusive owing to the difficulty of obtaining 
reliable data. The theoretical attack is similarly made 
difficult by the need to employ mixed Eulerian and 
Lagrangian coordinates. A discussion of electrostatic 
charging phenomena has also been included and it 
shows that in some systems it exerts an over-riding 
influence on the particle motion and momentum 
transfer. 


he various aspects of solids-gas flow which have been dis- 

cussed in the five previous articles of this series 
stricted to the behavior of single particles. Consideration will now 
be given to multiparticle systems in turbulent fluids in which the 
inertia of the fluid may be sufficiently large or the particles 
sufficiently small for the solids to have a turbulent motion of 
their own. This subject is one of great theoretical and analytical 
complexity and many of its aspects are still in a very preliminary 
stage of dev clopment. The experimental approach is equally 
difficult and precise methods of measurement of the physical 
properties of multiparticle systems are often lacking. ‘The 
discussion which follows reflects the uncertainties of the theoret- 
ical knowledge as well as the paucity of reliable data needed for 
a sound statistical characterization of the system. 


were re- 


Consideration will first be given to the mechanism of entrain- 
ment of solid particles in horizontal streams, and their effect on 
the fluid velocity profile. ‘This will be followed by a review of 
the few studies, both theoretical and experimental, which are 
concerned with the turbulent motion of a small particle in a 
turbulent fluid. No attempt will be made at this stage to deal 
with the momentum transfer aspects of dilute-phase pneumatic 
conveying, as this topic will be treated in a subsequent article. 
\ brief discussion of electrostatic charging phenomena has also 
been included because of its related interest. 
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X—MULTIPARTICLE BEHAVIOUR IN HORIZONTAL 
TURBULENT STREAMS 


Phenomena occurring in sediment-carrying horizontal streams 
should be of direct interest to workers dealing with dilute phase 
solids-gas flow. Vanoni stressed that although most of the 
mechanics of the transportation of sediment was investigated in 
connection with the control and dev elopment of natural streams, 
it is directly applicable to the flow of solids-gas mixtures in 
industry. The underlyi ing mechanisms of these two fields have 
much in common, but sediment transport theory has generally 
been ignored in investigations of solids-gas flow. 

When a bed of solids is submerged in a horizontal fluid 
stream of gradually increasing velocity, Font succeeding types 
of behavior may be discerned: (a) initially there is no solids 
motion; then (b) individual grains in isolated random areas of 
the bed begin a rolling or sliding intermittent motion in short 
steps; (c) some grains are actually lifted into the stream and 
they eventually settle out, colliding with particles which are 
thus dislodged from the bed and are in turn carried into the 
stream; this phenomenon which is referred to as saltation results 
in the transport of the material along the bed mainly in the form 
of continuously changing ripples so that saltation is best described 
by statistical methods; (d) in the final st age the grains are lifted 
and suspended by the cross-component of the fluid turbulence 
their trajectaries incre asing in length as the fluid velocity in- 
creases. The critical fluid shear stress on the bed surface that 
will just start particle motion is designated by the symbol rt... The 
solids movement in the horizontal solids- gas systems will exhibit 
a similar behavior if the particle concentration exceeds a char- 
acteristic value. 


Fluid Velocity Required to Initiate Motion 


White investigated the equilibrium of grains on the bed of 

a stream, and found that for spherical particles, tluid turbulence 

may reduce the required value of 7. by more than 40%. Rubey“ 

in an analysis similar to that of White, found the minimum drag 
required to be given by 

Rots = Cp ad Pp. U2/2 (1) 


where U’, is the fluid velocity at the grain and ts less than the 
average stream velocity. This expression Is sumilar to the 
general resistance equation in Ww hich Cp’ would correspond to 
the drag coethcient. This relationship was not successful tor 
beds made up of fine grains, however. 

Vanomi showed that since 7, is approximately proportional 
to U2 it would account for the rule of thumb long used by 
geologists that the mass of the largest particles that can be 
moved by a stream is proportional to the sixth power of the 
stream velocity. 

Shields” recognized that the fluid velocity at the grains was 
lower than the stream velocity, and he applied the Karman- 
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Nikuradse statement of the velocity distribution in a fluid stream 
to obtain the relationship: 


(p — po) gd = Q(U,d/v).. ; (2) 


Shields found that 7, tends towards a minimum when the particle 
Reynolds number, defined by (U;d/v) was equal to 10 while 
White® found that it was smaller for (U;d/v) > 3.5 than for 
(U,d/v) <2 


Bed Load Transport 


‘Bed Load” transportation theories have been presented to 
enable the calculation of the total sediment carried on or near 
the bottom of an horizontal stream. DuBoys proposed the 
expression: 

= Br.(tT — 7). (3) 
where w = mass flow rate of solids carried per unit width; 
B = sediment parameter. 
Values of B and 7, have been correlated and listed by Brown‘® 
but the applicability of the expression has been found to be 
limited. 


Einstein”, on the other hand, disputed this concept of a 
critical hydraulic quantity, and felt that the phenomenon is one 
that involves probability. He assumed the particles to migrate 
in jumps of length L proportional to the particle diameter. The 
number of particles migrating would be proportional to the 
number of grains in the surface layer and also to the probability 
that the local drag on the particles i is sufficient to lift them into 
the stream. 

Einstein defined two dimensionless parameters ®, and ®, 
which reflected the energy of particle and fluid movement 
respectively. These were given as: 

®, = (w/p) [p./(p — po)|''?/[gd3]'?............ (4) 
and 
en ae (5) 


and Einstein was able to show an experimental and analytical 
interrelationship between the two functions in support of his 
model. Einstein found that his relationship seemed to account 
for observations obtained from a system involving sand with a 
wide particle-size distribution. Vanoni, however, found it 
to be inaccurate when applied to systems of uniformly-sized 
particles. 

Kalinske®°'” proposed similar expressions, which gave fair 
agreement for systems of the latter type, of the form: 


= [(pd) VTo/po] | Dit</T.] ..-. (6) 
where ®, is a function which is related to the fluid turbulence. 


Solids Concentration Gradients in Horizontal 
Fluid Flow 


The solid phase of a horizontally-flowing solid-fluid system 


will often exhibit a concentration gradient which in the case of 


solids-gas flow may determine the flow characteristics and energy 
requirements of the system. An equation for the concentration 
of sedimenting particles i in a still fluid was proposed by Schmidt 
in 19252 and re-stated by O’Brien“ in 1933. It has the 
form: 


cu + D,(dc/dy) = 0 (7) 
where D, = a diffusion coefficient; 
y = the vertical distance from the lower boundary of 
the system; ; 
¢ = the concentration of solids; 
vw = the free-fall velocity of the particles. 


vy, would be more strictly defined as the fall velocity of the 
particles in a system whose turbulence parameters are equiy alent 
to those existing in the system to which the expression is being 
applied. This statement of dynamic equilibrium balances the 
settling rate of cv, with D, (dc, dy) the upward rate of movement 
of the solids due to the turbulent fluctuations of the fluid. D, 
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was found experimentally to be a function of y, and it wa 
measured by Kalinske®. 


Rouse“ derived an expression for the solids distribution 
in the following manner: 


By definition for a solids-free flow 
rt = p,D{dU/dy)..... (8) 
He then introduced two simplifying assumptions by saying that 
D, = D, and by assuming 7 to be a linear function of y, i.e. 


Cfte SA =P) cei - hem eee 

From the von Karman Law 
dU/dy = [Wr./pol/ky = U;/ky... Sit OS 
so that D, = ky U1 — (9/B)) occ ces (11) 


By substituting 7 into 11 and integrating, Rouse obtained: 





im a 
they a meet SE dace ee 
y h-wn 
where ¢, = concentration at level y,; and 
G1: = O/RU-)..55.. Ea F 8S) 


where %, is the settling velocity of the solids at level y;. 


Also, since To = poghS.... Sales ence 8 ae 
U, = (ghS)"....... a. bn ei 
and es BT isos Se hws eae oe (16) 


Vanoni“® and Ismail“? have experimentally verified the 
general form of this expression, but the experimental a, referred 
to aS (@1) exp, Was smaller than the theoretical value. In addition, 
the suspended matter was more evenly distributed than was 
predicted by the theoretical expression. 

Data subsequently obtained by Vanoni‘) showed agreement 
with Equation (13) when (a1) was adjusted to a proper value. 
Vanoni felt that (@i)exp, WaS not —- to the theoretical value 
because the assumption that D, D, did not take into con- 
sideration the smaller inertia of the fluid compared with that of 
the solids. He corrected for this by putting: 


Dy = a2Dy....... ; phenatoret aince th (17) 


:xperiment showed a» to vary from 1.0 to 1.5, and it seemed 
to be affected by the size of the particles. 

Ismail”) obtained the empirical relation D, = 1.5 Dy for 
0.10-mm. sand particles in water, and D, = 1.3 D, for 0.16-mm. 
particles all in the same system. Carstens“®) however doubted 
that a2 could be greater than unity, since D, is really the diffusion 
coefficient of the fluid. Since the density of the the fluid was at 
all times less than the density of the solids, D, should always 
be bigger than D,. Carstens cited evidence to this effect in his 
work and in that of Wagenschein“”. 

In reply to Carstens, Ismail claimed that the work of 
Sherwood and Woertz® and also that of Corcoran’? showed 
in effect that a2 could be greater than 1. Nevertheless the 
results discussed in Part III of this series, dealing with the 
resistance of a spherical particle in a fluctuating flow make it 
very difficult to uphold this argument. Ismail’s apparent ob- 
servations may suggest that the model employed is erroneous 
and the mixing length analogy cannot be applied to this system. 


The Effect of the Suspended Solids on the Fluid 
Friction Factor 

‘The effect of the suspended solids on the flow characteristics 
of the fluid was investigated by Vanoni“ and he found that the 
value of the Von Karman universal constant k was decreased in 
a fluid-solid system. In his experimental study, the flow rate 
and the depth of the channel were kept constant (i.e., the shear 
remained constant) and the solids content of the stream was 
increased. From Equations (8), (10) and (11), if 7, (and 
therefore U7) is kept constant, decreasing k would cause 
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(dU/dy) to increase so that Dy would decrease. Vanoni found 
that k would decrease to as low as 0.2 with an increase in 
stream loading. 

Vanoni noted that increasing (dU/dy) will tend to increase 
the mean velocity, and the friction factor for the channel will 
tend to decrease, a behavior which was confirmed by experi- 
mental observation. He attempted to explain the decrease in k 
by assuming that the presence of the suspended matter damps 
the turbulence since it must supply the energy per unit mass of 
fluid per unit time to prevent the solids from settling, and this 
would be proportional to the solids concentration and settling 
velocity. Results obtained with sedimenting systems have shown 
D, to increase with particle concentration, but Vanoni admitted 
that the process by which this occurred was not clear. In his 
observations the energy which would be required to suspend the 
solids amounted to only 3% of the amount required to overcome 
the friction in the channel, whereas the alterations to the friction 
factor were of a much larger magnitude. 

In a solids- liquid system the velocity of the particles with 
respect to the fluid, 1.e., their slip velocity, is very small relative 
to the values encountered in solids-gas systems. In systems 
with low slip velocities, the turbulence generated by this relative 
motion may be negligible in comparison to the turbulent energy 
being used up to suspend the solids phase, and the net effect of 
the solids is to exert a damping action on the fluctuating velocity 
of the conveying medium. In a solids-gas system, ‘the finite 
relative velocity could generate turbulence which will be dis- 
charged from the wake of the solids in much the same manner 
as in the production of turbulence by the placement of a wire 
mesh over the cross sectional area of the moving stream. It 
may very well be that in this case the turbulence generated by 
the particles would exceed that which is withdrawn from the 
system in the process of the suspension of the solid phase. Other 
aspects of the problem of the reaction of the fluid phase to the 
presence of the solids will be taken up in discussions which 
follow. 


XI—ALTERATION OF FLUID VELOCITY PROFILES 
DUE TO PRESENCE OF PARTICLES 


Laminar Flow 

The presence of solids particles in a laminar stream flowing 
through a pipe will cause an alteration in the characteristic 
parabolic velocity distribution. Brinkman? deduced an expres- 
sion for the flow pattern in a vertical column by analogy with 
Poiseuille’s law to give: 


Pw dp Tol(y/F ” a a 1] 
Ue —— = pe - (18) 
M dx L I, a/VP m 


where J, = Bessel function of zero order and first kind and of 
imaginary argument; 
y = distance from the axis; 
and P,, is the permeability of the system defined by: 


Pn = (77/18) {(3 + 4/b. — 3 [(8/ho) — 3]2}....... (19) 


Plots of the velocity profile for different values of 4, showed 
that the distribution becomes less clongated as @, is increased 
so that it is almost uniform at P,, = 1072 ft.2 This result has 
not as yet been confirmed by experiment. 

If the particles become sufficiently small, the multiparticle 
system will begin to assume the properties of a non-Newtonian 
fluid, and the results of the investigations of these systems is 
not of direct interest to the solids-gas flow problem. 


Turbulent Flow 

The results obtained with sedimenting streams give some 
indication of the alterations in the velocity profile of the stream 
caused by the presence of the solids. A better analysis of 


Vanoni’s data than had originally been given by Vanoni himself 


was suggested by Hunt) who showed the velocity distribution 
in a particle-laden stream to be given by: 
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Umax ze U 
(ghS)'2 k ke 


1 oc 1/2 
= —~4(1—y/h)! 2 + kaln se 





.(20) 


where ky = constant of integration: 


Umax = Maximum velocity of the stream: 
y = distance along the y axis, i.e., the axis transverse to 
the flow; 
h = total depth of flow. 


[he particle concentration effect acts by altering the value of 
k, and this must be determined by experiment. 

Laursen and Lin‘ felt that bed conditions alone govern the 
variation of ks which is used in defining the velocity profile as: 


SS oie oleae Sod .. (21) 


They believed that the suspended sediment has no effect on the 
production or distribution of turbulence, and that the phenomenon 
is similar to that noted by Nikuradse in connection with the 
change of the velocity distribution in a pipe with the roughening 
of the pipe wall. 

In a communication to Laursen and Lin, Ismail“? noted that 
their conclusion was based on values of the function &3; and on the 
contention that as the velocity increases, ks tends to the value for 
solids-free fluid. Ismail reported detailed experimental ob- 
servations which suggested that this was not correct. 

Very little is quantitatively known about the effects of the 
solids phase on the velocity profile of a gas flow other than that 
the pressure drop in two- phase flow which is contributed by the 
gas phase cannot be determined from the conventional friction 
factor plots for solids-free flow. No method has as yet been 
presented for the calculation of this quantity, and this has 
hindered the development of a general theory for solids-gas 
flow. Coulson and Richardson» pointed out that small con- 
centrations of solids in solids-liquid pipe flow have caused a 
reduction in the overall pressure drop encountered. Experimental 
measurement of the effects of model fibres in the flow character- 
istics of fibre suspensions carried out in this laboratory®® have 
shown that for short fibres (less than one millimeter in length) 
no flocculation occurred and the reduction in pressure drop 
(below the value for water alone) increased as the fibre con- 
centration increased. For longer fibres, flocculation did occur, 
which was accompanied by a marked increase in pressure drop 
over the value for water alone, owing to the establishment of a 
laminar regime of flow. Turbulent conditions at higher velocities 
caused the disruption of the floes, and the pressure drop fell once 
more below the water value. 

Knowledge of the velocity profile of a solids-fluid system 
would allow the measurement of the pressure drop due to the 
fluid phase. However, velocity profile data for solids-gas systems 
are lacking, mainly because of the difficulty in employi ing a 
pitot tube in two- phase flow. Dussard and Shapiro®? have 
recently designed a modified form of the pitot tube for the 
measurement of the velocity profiles of particle-laden streams, 
and this instrument may become an important tool for the 
investigation of solids-gas systems. 

Zenz*), in a study of solids-gas flow in pipes found that in 
the case of one material there was no appreciable difference in 
the pressure drop for the flow of the solids-gas mixtures and the 
flow of solids-free air. This suggests that the presence of the 
solids altered the turbulence characteristics of the fluid so as to 
lower the frictional retardation which it encountered at the pipe 
wall. 

Segler® recorded the velocity distribution in a 9-in. pneu- 
matic grain conveying tube through which varying concentrations 
of wheat grains were being transmitted. The velocity profiles, 
shown in Fi igure 1, were found to vary with the solids concen- 
tration and the greatest change occurred in the vicinity of the 
pipe wall. Increased loading ‘elongated the profile and decreased 
the velocity gradient at the wall. No attempt was made to 

calculate the pressure drop due to the movement of the gas 
phase with an increase in solids concentration. This observation, 
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Figure 1—Velocity profiles of air stream containing solids. 


if correct, would be of great assistance in interpreting and 
explaining data obtained in investigations of dilute phase solids- 

as flow, since it suggests that the ‘fluid pressure drop is at times 
almost half of what would be obtained with solids-free flow. 


XII—THE TURBULENT MOTION OF A SMALL 
PARTICLE IN A TURBULENT FLUID 

The previous discussion has shown that the alterations in 
the mean and turbulent characteristics of a flow by entrained 
solids will be related to the motion of the constituent solid 
particles. Three particle regimes occur in solids-gas systems. 
If the velocity of the particle is such that its mean relative 
Reynolds number is in excess of say 100, its motion will be 
essentially linear, and its effect on the free-stream turbulence 
will be similar to the action of a screen placed normal to the 
flow. The investigations of Schwarz‘° and those of the present 
authors‘, discussed in Part V of this series, were performed in 
this region. At the lower Reynolds numbers of this regime, 
the wake turbulence contributed to the flow will be associated 
with high wave numbers, and it will hasten the dissipation of 
the fluid turbulence. At higher Reynolds numbers, the wave 
numbers will become sufficiently low so that there will be a 
genera! increase in the stream turbulence level. Even when the 
gross level of turbulent energy is reduced, the vorticity shed by 
a particle would have frequencies which would be the most 
effective in accelerating the spectral momentum transfer in the 
wake of the neighbouring downstream particles. 

A second particle regime occurs when the particle Reynolds 
number is less than unity. Here the particle mean velocity is 
essentially equal to that of the free-stream, and the particle 
acquires a turbulent motion which is a response to the fluid 
turbulence in its immediate vicinity. Intuitively it would be 
expected that as the size of the particle decreases or its density 
approaches that of the fluid, the particle motion will follow more 
faithfully the motion of the fluid which it displaces. The discus- 
sion in this section will deal with this second regime. 

Separating regime | and 2 will be a third intermediate state 
of motion where the particle has both a mean velocity with 
respect to the fluid, i.e., a finite mean Reynolds number, and a 
turbulent motion which is a response to the ambient fluid turbu- 
lence. The theoretical description of this regime will have to 
await further mathematical progress with mixed Lagrangian and 
Eulerian statistics. 


The simplest model from which the problem of the turbulent 
motion of a particle can be broached is a stationary one in 
which the particle is suspended in an infinite fluid element 
which is oscillating with a periodic sinusoidal motion, and a 
description is sought of the particle fluctuation frequency and 
amplitude. The theoretical and experimental treatment of this 
system has been discussed in Part III of this series”. This 
stationary model has been adapted to the turbulent fluid case by 
giving the coordinates of the stationary system the velocity of 
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the particle being considered. A solution of this type was 
obtained by Tchen‘*?) who assumed the particle to be in a 
periodic field defined by: 


A similar solution was obtained by Liu“® who used a 
generalized harmonic approach in which the velocity of the 
particle i is represented by a Fourier integral. This allowed him 
to obtain an expression of the function ~(w) which denotes the 
ratio of the sinusoidal component of the velocity of the spherical 
particle to the corresponding component velocity of the fluid 
clement at the position of its centre, when the spherical particle 
The function is 





2 


is absent. In other words, ¥(w) = u’,2/u'?. 
simplified for the following special cases: 


(1) Non-viscous fluid: If the viscous effect is negligible com- 
pared to the inertia effect: 


BRB eco ee ca Hae he sealer HA (23) 
(2) If (rw!!2)/(v'2)<K 1 
Y(w) = 1 — (1)/[1 — (9/2) (ipsv)/(p — po)wr?))..... (24) 


and if (wr?)/(v) >>(p.)/(p — po) 


¥(w) becomes very small and the particles remain nearly at 
rest as the turbulent flow beats on them. 


(3) As the sphere radius r decreases to a very small value, 
wW(w) tends to 1 and this will occur when: 


rF#<K (9/2) [(p.)/(p — po) (v/w)].. 2.0... (25) 


Solutions of this type require that the particle diameter be 
much smaller than the Eulerian microscale, and that the particle 
always remain associated with a given packet of fluctuating 
fluid. Even then it is doubtful whether one can assume that 
Stokes’ Law applies since a fluid element in a turbulent field 
will be undergoing continual distortion. Corrsin and Lumley“ 
have shown that a further weakness in the derivation develops 
when it is applied to a turbulent fluid. If the coordinates are 
fixed to the particle, the inhomogeneity of the particle movement 
will be related to time alone, but the fluid velocity inhomogeneity 
will still remain a function of space as well as time. Lumley“ 
has elaborated on the extreme mathematical difficulties which 
occur when this factor is incorporated into an attempted solution. 


Friedlander“ has used a mathematical approach similar to 
that used in the calculation of a nian movement to obtain a 


relationship between u’,? and w’,?. Here again it was necessary 
to assume that Ke — 0, and a the particle diameter was 
much smaller than the fluid Eulerian microscale. In addition, 


(Nu? , p./m) was restricted to being less than 1, and the 
Stokesian drag was assumed to be valid for the eelens resist- 
ance. For very small particles entrained in a liquid, Friedlander 
assumed that the Lagrangian velocity correlation would suffice 
(i.e., the particle would remain associated with a given packet 
of fluid so that correlations as a function of time — could be 


employed) and an approximate expression for u’,? = u’,* — u’? 
was found to be given by: 


u'? = 2(1 — p,/p)? u’P/AB?..... fade 
where A is the Lagrangian microscale, in sec. 
and 8 = 3apd/m......... sin 0 SORE 


For the case of solid particles of unit specific gravity in a gas 
flow, Friedlander applied the simultaneous experimental meas- 
urements of the Eulerian and Lagrangian correlations given by 
Mickelsen‘*”) for an eight-inch duct which were discussed in 


Part I. This led to an approximate expression: 
40’ 2/ulf m s/(8 +B)... - 6c ccc cease (28) 
where 2 is an experimental constant. A plot of w’,?/u’? as a 
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the former to rise from zero at d = 0 to 0.2 atd = 


, 


The ratio approached 1, i.e., “’)? — 0 as d approached 100 
microns. The solution was not suggested for use above 10 
microns however, since there would be a mean relative velocity 
between the particle and the fluid, and the particle would have 
a finite mean Reynolds number. In addition, a plot of Nu’, the 
Nusselt number for mass transfer, was shown simultaneously 
against d. Nu’ rose from a value of 2 at d = 
about 4.2 at d 100 microns. Friedlander commented con- 
cerning the latter finding that the assumption sometimes made 
of Nu = 2 (stagnant case) for particles carried along by a 
turbulent fluid may be in error. Recent studies by Dlouhy (38,39) 
on the other hand, have shown that small droplets (from 10 
to 40 microns) will evaporate in turbulent streams with Nu = 2. 
It should be noted, however, that the Mickelsen’s data, on which 
Friedlander based his calculations, may be unreliable“. 

By using similar simplifying arguments and restrictions and 
by neglecting the effects of acceleration Soo” obtained solutions 
of the particle turbulence parameters as a function of those of 
the fluid. The solutions were given as a family of curves showing 
the ratios of the intensities, scales and diffusivities of the particles 
to those of the fluid as a function of: 


© Re d p (29 
Kreis RNase 29) 
where Re = dvu? Deh Mesa ridin weeta ets (30) 


and is the Lagrangian microscale of turbulence. 

The relationships which he obtained were not applicable to real 
systems, and Soo decided that because of the difficulties involved 
in obtaining adequate solutions the interrelationships would have 
to be obtained experimentally. 


In what is perhaps the most sophisticated investigation of 


solids-gas flow to date, Soo developed experimental and com- 
putational methods which allowed him to measure the statisti- 
cal properties of both the particulate and entraining fluid 
motion“?:49, The experimental system consisted of a 3-in. 
square duct with a 10-ft. horizontal plexiglass test section. 
Glass beads, 105 to 210 microns in diameter, were fed at rates 
up to 0.5 Ib./min. into the air stream which was varied up to 
85 ft./sec. A thin vertical sheet of light was used to illuminate 
the system and the particle motion was recorded by stroboscopic 
photography. By drilling equal-sized holes through the sequence 
of particle images on the photographic plates, Soo used an optical 
autocorrelation technique to measure the Lagrangian scale and 
intensity of the particle motion. The average Lagrangian 
properties in the central core of the fluid phase were estimated 
with the particles present by using a helium diffusion technique. 

The particle Reynolds number was in all cases less than 

), and there was no significant alteration of the average fluid 
eas properties by the presence of the solids. Two 
remarkable results were obtained. In the first place, although 
the fluid turbulence would be expected to be roughly isotropic 
in the central third of the duct in which the measurements were 
taken, the particle intensity in the y direction No’, was con- 
siderably smaller than the particle intensity Nu’)? in the axial 
direction. Soo feels that this is a direct result of the energy 
required to overcome the gravitational forces tending to cause 
the particles to settle out. This could have been confirmed by 
using a horizontal sheet of light and measuring the z component 
of the particle intensity. 

The comparison of the particle intensities with the fluid 
intensities provided an even more startling result in that the 
former were in general considerably greater than the stream 
turbulence. The theoretical solutions obtained by Soo had 
shown the ratio: 


T = (uP — u'p)/(u'P)..... .»» (31) 
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function of the particle diameter for a mean air velocity of 


0 to a value of 








to increase from zero at ®, < 0.1 and then increase to 1 at 
®, > 50. The experimental results all showed negative values 
for this ratio, at times going as low as —12. Using a more 
valid set of starting equations, than in his original solution, 
Soo attempted to obtain an approximate solution for the particle 
turbulence parameters for varying particle distance from the 
duct wall), This showed that the effect of the wall was to 
displace the T versus ®, curves upwards rather than downwards, 
and at no time were there any negative values of T predicted. 

Soo then attempted to explain the results as being caused by 
the increase in fluid turbulence, with increasing distance from 
the duct centre, discussed in Part 1%. If a particle was in an 
off-centre region where the ratio of the fluid turbulence at that 
point to the ‘turbulence at the duct centre was say 2, and if T 
was actually zero, the value of T which would be calculated using 
the central fluid intensity would be —3. Soo’s observations were 
taken in the central third of the pipe so that it is difficult to see 
how the higher particle intensities could be transported into this 
region by particles which were in the vicinity of the wall. Even 
if this argument is granted, however, it still does not explain 
the large particle i intensities which were observed at the centre. 
For example, the experimental data at a mean fluid Reynolds 
number of 40,000 and ®, = 10 show T values of —11 whereas 
the largest negative value of T which could be attributed to the 
maximum intensity in the whole duct cross-section (i.e., the 
sharp rise near the wall) would only give a value of —1. 

The reason for the particle intensities being higher than the 
fluid is difficult to explain, since Soo’s careful experimental 
techniques leave little doubt that this behavior actually occurs. 
It must be that the particle velocity fluctuations are a response 
to forces other than the fluid turbulence. The generation of 
static charges, which is to be discussed in the following section, 
is very difficult to avoid in solids-gas systems, even when 
elaborate grounding and humidifying procedures are employed, 
and mutual particle repulsion would perhaps contribute to the 
fluctuating motion. A second possibility could be that the drag 
forces on the particle in a turbulent packet of fluid are more 
complicated than those given by the Stokesian solution employed 
and that they themselves may act in an unsymmetrical and 
unsteady manner. In spite of the difficulty of offering a plausible 
explanation for the results which he obtained, Soo’s researches 
are representative of the type of investigations which will 
eventually provide a scientific understanding of dilute-phase flow. 

Peskin“**-45.46.47) has raised the possibility that fluctuating 
pressure forces are exerted on each particle in a solids-fluid 
system by its neighbors, similar to those predicted by potential 
flow theory. Konig(*®) for example has shown that two > particles 


moving in a ficticious inviscid fluid will exert a force F on each 
other given by: 

= (3/2) (p,/m) Ur? (r8/I|5)0. ee (32) 
where / is the displacement vector between the particles. 
Peskin assumed that this force term can be extended to non- 
steady motion and can be applied to particles which are in 
turbulent motion in a turbulent velocity field. In order to 
proceed, he assumed that the solids will not affect the fluid 
turbulence, and that the fluid drag on the particles is Stokesian. 
The extremely difficult task then becomes to evaluate the total 
interparticle pressure force acting on a particle due to its neigh- 
bors. If the assumptions made are granted, the force is given by: 


N 
Ft) = -* OP Fe ms wastes (33) 
2m gag |4ii° 
where w’ is the intensity of fluid motion relative to the particles. 
An important consequence of the solution of (33) is that there 
is a particle diffusivity due to the random pressure forces, but 
that these are felt only by slow particles, and can be ignored 
as the particle velocity is increased. Heavy particles will, of 
course, tend to have very low diffusivities. 


117 









The method of dealing with random interparticle forces 
presented by Peskin should be directly applicable to the case 
of electrostatic repulsion in charged systems. 


The Effect of the Solids on the Fluid Turbulence 
Level 


Measurements cannot be taken with a hot-wire anemometer 
in a particle-laden stream so that indirect turbulence measure- 
ments must be obtained from diffusion experiments. The 
information which is available from these results is interesting, 
but they are limited as to detail since they provide intensity and 
scale values which are the av erage across the diffusion breadth. 
Particles with sufficiently high Reynolds numbers would be 
expected to contribute turbulent energies in a wave number 
range which may have a maximum effect on the downstream 
particles, and yet their contribution to the overall level of the 
stream could be insignificant. Spectral density measurements 
would in addition provide an insight into the energy consumed 
in suspending the particles but the development of instruments 
to obtain this information will be extremely difficult. 

The tracer diffusion technique has been employed by 
Hanratty, Latinen and Wilhelm“® to obtain turbulence meas- 
urements in fluidized beds. Their results with water-fluidized 
3-mm. glass spheres show that the intensity decreased steadily 
with increasing voidage, and the intensity values at a voidage 
of 90% were ‘essentially equal to the particle- -free case. The 
scale value from 70% voids upwards also decreased steadily to 
the empty tube result. 

Kada and Hanratty“*® have recently extended the use of the 
tracer technique to study the effect of solids on the fluid turbu- 
lence of cocurrent solids- liquid systems. Here the difference 
between the solids and fluid velocities is very small and would 
be extremely difficult to measure. The authors estimated it by 
assuming that the free-stream turbulence would not affect the 
particle boundary layer and wake structure, and the particle 
Reynolds number range calculated by this means was 0.66 to 
18.3. The product of the square of the relative intensity and 
the particle Reynolds number was in all cases much less than 
45 so that the turbulence effect should be small, At the lower 
particle Reynolds numbers, the solids did not affect the fluid 
turbulent diffusion coefficient up to a solids concentration of 
1%. With a particle Reynolds number of 18.1 however, there 
was a 20% increase in the oe coefhcient with a solids 
concentration of only 0.15%. Up to the limiting concentration 
of this study (2.5% solids) the turbulent diffusion coefficient 
increased sharply for low fluid Reynolds numbers, but increased 
very gradually with a Reynolds number of 50,000. The accumu- 
lation of data of this type is very important for an understanding 
of solids-fluid systems, particularly since the results are not 
complicated by electrostatic effects. 

The only results available in dilute-phase solids-gas flow are 
the ones discussed previously given by Soo‘), Unfortunately 
these were taken at particle Reynolds numbers less than 10, 
where there would not be any contribution from the wake, 
nevertheless they indicate that any damping action if present 
was restricted to the low energy, high frequency end of the 
spectrum. 


XIII—ELECTROSTATIC CONSIDERATIONS 


Solids-gas systems differ from their solids-liquid counterparts 

in that they tend to generate static electricity. This effect may 

‘ play a significant role in obscuring data by introducing force 
terms which are not accounted for in the theoretical equations 
being tested. 

Coulson and Richardson‘*®) noted that electrostatic charge 
accumulation in a solids-gas system may increase the pressure 
drop by a factor of 10. They feel that considerable further 
research is needed for a better understanding of the phenomenon. 
Lewis, Gilliland and Bauer‘) found that fluidized particles 
which were treated so as to minimize electrostatic charging 
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gave characteristic curves which differed in shape from those 
obtained with untreated particles. 

Mehta‘? studied the flow of glass spheres suspended in an 
air stream moving through pyrex pipe. He observed extensive 
charging between the particles and the air, and the particles 


and the pipe wall. Long axial arcs were formed every few 
minutes and these were succeeded by pressure waves which 
caused error in the manometer readings. The discharge caused 
the combustion of traces of oil entrained in the air, depositing 
carbon on the pipe walls. Mehta noticed that the pattern of 
flow, i.e., the distribution of the solids phase oven] the flow 
cross-section could be altered at will by altering the | ppsition of 
the grounding wires, and no amount of grounding wis able to 
mitigate the charging action, 


Mitlin®*®) referred to the ‘‘anomalous behavior’’ of sand 
particles in a metal pipe. ‘The pressure drops encountered 
would gradually increase with time to a constant value. It 
reverted to its initial value after the flow had been discontinued 
for a lengthy period of time. Mitlin felt that this could be 
explained by the electrostatic charge residing on a very fine 
dust layer on the inside surface of the tube. This phenomenon 
is held to account in part for the increase of pressure drop with 
solids concentration. In their study of the elutriation of fine 
particles from dense-phase fluidized beds, Chin-Yung Wen and 
Hashinger“® recently found that reproducible data could be 
obtained after “the electrostatic charge was sufficiently elimi- 
nated from the column” by wrapping grounding wires around 
it and also introducing a grounding wire into the bed. They 
questioned the validity of a previous study®® on the basis of 
the possible existence of electrostatic charges. 

Culgan®® found that in a short time the metal wall of a 
conveying system carrying plastic particles will become coated 
with a hard lay er of this material. This could result in a system 
similar to the one noted by Mitlin, and it would be difficult to 
define the electrostatic generation components without knowl- 
edge of the actual surface composition of the pipe. 

Richardson and McLeman‘*”) have recently reported a very 
detailed investigation of the development of electrostatic forces 
in solids-gas systems which reflects an increased awareness of 
the importance of this phenomenon. The solids were conveyed 
in a horizontal 1-in. diameter brass pipe which was 114-ft. long. 
When 0.05-in. sand particles were recirculated through the test 
loop, the relative velocity between the solids and the conveying 
gas gradually increased to an equilibrium value. The charge 
buildup was not associated with the particles since the substitu- 
tion of fresh sand had no effect. Oddly enough cleaning the 
pipe surface, grounding or charging it produced no visible effect. 
Correlations were obtained between the pressure drop and the 
relative velocity increase, but no attempt was made to relate 
their variations with the time of conveying. Surprisingly, no 
correlation was found between the charges on the particles and 
pipe surface on the one hand and the relative velocity and pressure 
drop on the other. Richardson and McLeman ‘feel that this 
supports the existence of a charged dust layer at the tube wall. 
With the sand system, the dust required an Opposite charge to 
the particles, and the resulting attraction for the wall region 
retarded the particle motion. When Perspex particles were 
used, the pressure drop decreased with time, and in this situation 
the charged dust layer was felt to have acquired a similar charge 
to the particles w hich caused them to move near the axis of the 
duct. 

Khudiakow and Chukhanov“**) found that electrostatic re- 
tardation of particle movement could become large enough to 
completely halt the flow of the solids phase in a vertical solids-gas 
system. High speed photographs had indicated that lesser 
amounts of charging had caused the particle motion to be 
“rather drastically changed”’ 

Electrostatic discharge is sometimes a hazard which prevents 
materials such as sulphur from being conveyed pneumatically. 
Sinha” as well as Gasterstadt® reported the generation of 
large sparks in industrial solids-gas systems. 
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There has been no attempt at a quantitative evaluation of 


the effect and extent of electrostatic charging occurring with 
motion of particles in a pipe. 

Green“ developed an equation for the motion of a charged 
smoke particle introduced into a stream of gas moving with 
laminar flow between two vertical plates maintained at a potential 
gradient E,: 

v = velocity of movement of the particle to one of the plates 


ne i) . (34) 


This indicates the type of expression which may have to be 
incorporated into the momentum equations defining solids-gas 
systems. 

It would be advantageous to be able to predict the electro- 
static generation qualities of the various combinations of solids- 
phase and pipe wall materials encountered. This would however, 
necessarily involve a better understanding than is presently 
available of the mechanism by which friction generates static 
electricity. 

Thomson *) held that frictional electricity is really merely 

contact electricity, the contact being made good by considerable 

pressure. Vieweg“** cited the Helmholtz theory of frictional 
electricity which postulated that the fundamental cause is a 
contact difference of potential which exists between any two 
materials. When two uwnlike surfaces are brought together, 
there is a concentration of charge at the surface of contact. 
On separation, if either is a nonconductor, the contact potential 
is greatly magnified and both surfaces are charged. If both 
surfaces are conducting, the charges are said to neutralize on 
separation in the last instant of charging. 

A more recent paper by Harper” showed that it is erroneous 
to think that no rubbing potential will occur if both components 
of the system are conductors. The charging of metal powders by 
blowing against metal surfaces is found to be comparable to that 
obtained when insulating powders are used. Experiments indi- 
cated that one cause for the electrification of metal metal 
surfaces is a manifestation of their contact potential, and Harper 
felt that these must always be present. The magnitude of the 
charge generated may be calculated from the contact potential and 
the topography of the surface involved by the use of a quantum- 
mechanical theory based on the original Volta-Helmholz hypo- 
thesis. Subsequent experiments by Harper‘® seem to show 
that for light contact with metals, hydrophylic surfaces such as 
glass acquire a considerable charge, whereas hydrophobic sur- 
faces such as polythene, amber, polystyrene, nylon, Perspex, 
silicone and Teflon do not. 

The mechanisms of the frictional electrification of particles 
on dispersion into a solids-gas system are still little understood ® 
and the data obtained are not reproducible and are often incon- 
sistent and even contradictory. 

The charges on sprays of electrified particles were studied 
by Sachsse“*” and Loeb‘*® for systems of dust particles, and by 
( hapman‘*®) for atomized liquid droplets. These indicate positive 
and negative charges to be present in equal numbers. In some 
cases the charge was proportional to the surface area, whereas 
in others it was found to vary as the particle radius. 


Kunkel“® observed hotographically the trajectory of 
P g j 


charged solid particles in an air stream. It was found that 
if a mass of particles were separated to form a cloud, practically 
all the particles would become charged, but the cloud as a whole 
would be neutral. If the cloud would come into contact with a 
surface, an asymmetry of charge might or might not occur. 
Asymmetry was usually observed in systems of metals i 

combination with insulators where the insulator usually receiv ie 
the negative charge. Kunkel found that in some cases, i.c., talc 
particles on a metal, a highly charged dust coating would cling 
to the metal. A surprising conclusion of the trials was that 
humidity does not affect separation charging. No change was 


observed as the relative humidity was increased to a value of 


90% at a room temperature of 20°C. Also, fluid turbulence was 
found to have little effect on the charging of the particles. 
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CONCLUSIONS 


Progress is being made with attempts to predict theoretically 
the turbulent behavior of solid particles as a function of the 
fluid turbulence. These researches have uncovered such diffi- 
culties that experimental observations will have to be relied on 
for any useable information. Ingenious empirical methods are 
being dev eloped and preliminary results suggest that the starting 
equations used in the theoretical derivations may not take full 
account of the various forces acting on the particle. 
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Nomenclature 
Any consistent set of units may be employed. Those given 
are merely illustrative. 


a = tube radius, ft. 

a’ = hydraulic radius of stream, ft. 

A = amplitude, ft 

B = sediment parameter, ft.? sec.3/Ib. mass. 


solids concentration, Ib./ft.8 


ll 


Cp = drag coefficient, dimensionless. 

Cp’ = coefficient of minimum drag to initiate motion, dimen- 
sionless. 

d = particle diameter, ft. 

D; = diffusion coefficient for fluid, ft.2/sec. 

D, = diffusion coefficient for solids, ft.2/sec. 

e = electrical charge, coulombs. 

E, = potential gradient, volts/ft. 

E = energy gradient in direction of flow, ft.-lb./ft. 

F,,F, = experimental factors defined in text, dimensionless. 

F = interparticle force between two moving particles, 
poundals 

h = total depth of flow, ft. 

g = local acceleration of gravity, ft./sec.? 

g. = dimensional constant, (ft.) (Ib. mass)/(sec.)? (Ib. force) 

Re = Bessel function, Equation (18) 

= von Karman constant, dimensionless 

k,k2,k3, etc. = experimental constants, dimensionless 

L = length, ft. 

l = interparticle distance, ft. 

m = particle mass, lb. 

Nu = Nusselt number, dimensionless 

Nu’ = Nusselt number for mass transfer, dimensionless 

ne = number of electrical charges, dimensionless 

p = pressure, lb. force/ft.? 

Pu = permeability of system, ft.? 

r = sphere radius, ft. 

R = fluid drag force on particle (lb. mass)/(ft./sec.?) 

Re = particle Reynolds number based on mean motion, 
dimensionless 

Re = turbulent Reynolds number (dVu’? p, u), dimensionless 

S = slope of stream bottom, dimensionless 

x = ratio defined by Equation (31) 

uy = fluid element velocity, ft./sec. 

u's = fluid fluctuating velocity, ft./sec. 

u's,v’» = particle fluctuating velocities, ft./sec. 

u', = relative particle fluctuating velocity, ft./sec. 

U = mean fluid velocity relative to particle, ft./sec. 

U, = critical fluid velocity required to initiate particle motion, 
ft./sec 

U, = mean fluid velocity, ft./sec. 

U; = friction velocity = (7,/p,)', ft./sec. 

v = particle velocity, ft./sec. 

vy = free-fall velocity of solids, ft./sec. 


mass flow rate of solids carried per unit width, (Ib. 
mass )/(sec.) (ft.) 

distance along axis transverse to flow, ft. 

experimental constant, sec.~! 

@1,Q2,Q3, etc. = experimental factors defined in text 


ll 


4 
Il 


B = 3rpud/m, sec.—! 

0 = time, sec. 

dA = microscale of turbulence, ft. or sec. 
im = fluid viscosity, lb./(sec.) (ft.) 

v = kinematic viscosity, ft.2/sec. 

Po = fluid derisity, Ib./ft.* 

p = particle density, lb. /ft.' 

iy = bulk density, lb./ft.* 
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T = fluid shearing stress in solids-free flow, (lb. mass) 
(ft./sec.?)/(ft.?) 
T. = fluid shearing stress required to initiate particle motion, 
(lb. mass) (ft./sec.?)/(ft.?) 
Te = mes shearing stress at bed surface, (Ib. mass) (ft./sec.?)/ 
(1t.*) 
@; = fluid movement parameter, Equation (5), dimensionless 
@, = Kalinske function, Equation (6), dimensionless 
@, = particle movement parameter, Equation (4), dimen- 
sionless. 
®, = Soo parameter, Equation (29), dimensionless 
d» = volume concentration of solids, ft. of solids per ft.® 
¥(w) = ratio of u’,? to u’?, dimensionless 
@ = angular frequency, radians/sec. 
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Optimal Bypass Rates for Sequences 
of Stirred Tank Reactors’ 





RUTHERFORD ARIS? 


A general method of finding the optimal volumes, 
temperatures and bypass rates for a sequence of 
stirred tank reactors is developed. After a discussion 
of the method of dynamic programming, the example 
of a sequence of adiabatic reactors of equal volume 
with a single reaction is used as an illustration. A 
simple graphical construction for determining the 
optimal bypass rates is given for this case. The limit- 
ing cases of the optimal control of a batch reaction 
by addition of cold feed is also solved. 


W: propose to consider several design problems for se- 
quences of stirred tank reactors which arise when we 
require the operation to be optimal. It is well known that for 
certain reactions a higher yield may be attained by using two 
reactors at slightly different temperatures"? 3), Denbigh has 
given a very ‘striking illustration of this for a system of four 
reactions and while so dramatic an increase in yield (25% with 
one reactor, 57% with two) is rare, even small increases in 
efficiency can represent considerable gains in profitability. The 
problems we wish to take up here are concerned with optimal 
bypass rates; other problems with stirred tanks have been 
dealt with elsewhere”. 

The general method to be used is that of dynamic pro- 
gramming which may briefly be described as follows. Consider 
a sequence of R stirred tanks (shown schematically in Figure 1) 
which we will number from the end to the beginning, so that 
the feed enters reactor R, passes to reactor R—1 and so on, 
leaving as product from reactor |. At each stage we have 
certain design or control variables such as the temperature, 
pressure or holding time, which must be chosen in some sense 
optimally. The choice of design variables at stage r affects the 
course of the reaction in that reactor, and so the transformation 
of the state of the process stream that this reactor effects. The 
state of the stream may be specified by the concentrations, 
temperature ete. and we denote the set of values of these quan- 
tities in reactor r and in the stream leaving it by p,. Thus the 
state of the feed is denoted pg+: and that of the product by pi. 
If the design variables at stage r are denoted by q, the equations 
governing the behavior of the reactor may be written 


pr = S(pritj Gr)... (1) 


This is perhaps a little abstract but the definite equations will 
be given later and for the moment this formula only means that 
when we have chosen the design variables q, we know what 
the — state of the reactor r, p,, will be for any feed state, 
p11. This knowledge may be through a thorough understanding 


1Manuscript received November 1, 1960; accepted February 27, 1961. 
2Department of Chemical Engineering, University of Minnesota, Minne- 
apolis, Minn., U.S.A. 
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of the kinetics of the reaction and behavior of the stirred tank, 
or it may be purely empirical, a summary of operating experience. 

Let us define a profit function P(p; is- qx) which ex- 
presses our estimate of the gain acruing from the reaction. 
In its simplest form P might simply be the yield of one particular 
species, but at a more realistic level it would take into account 
the cost associated with the design variables, g,. In any case, 
given the feed state pe+1, the profit P can be evaluated by means 
of Equation (1) as soon as the design variables qi,...qg have 
been specified. We want to find the ‘optimal design policy, that 
is the set of values qi,...qx that make P as large as possible. 
There will normally be restrictions on these design variables, 
as for example when the temperature must be beneath some 
upper limit, and a design which satisfies these restrictions is 
called admissible. When the admissible policy that maximizes 
P has been found, this maximum value is a function only of the 
feed state and number of stages, so that we may write 


R( Pri) = Max P...... (2) 


We must also emphasise that the optimal policy is a function 
of R and Pr+i and this we do by referring to it, where necessary, 
as the optimal R-stage policy with respect to the feed state pr+i. 

We will consider a special form of the profit function that 
suffices for all our ee namely 


P = V(pi) — Vibr+i) — Clq) —..-C(qr) 
R 
= > {V(p-) — Vi peri) — Clq)}...-. ipo 
1 


Here V(p,) represents the value of the process stream in the 
state p, and C(q,) represents the cost of operating with con- 
ditions g,. Each term of the sum in Equation (3) represents 
the net profit from one reactor as it increases the value of the 
stream from V (pr+1) to V (p,) at a cost C(q,). 
Let us build up the optimal R-stage policy a stage at a time. 
If R = 1 then 
fi(ps) = Max {| V(pi1) — Vip2) — C(q)},. (4) 


where this maximum has to be obtained by a right choice of qu 
and p; is given by 
bp: = SC ps3 q)...-.. (5) 
This problem presents little difficulty for the number of design 
variables in qi is quite small. Now let this one stage be the 
second of two stages from which the profit will be 


{ V(p2) — Vi bs) — C(q2)} + | Vp) — Vi b2) — C(q)}. (6) 





Figure 1—The general stage-wise process. 








But whatever choice of q2 is made we certainly cannot do better 
than to make the optimal choice of g; with respect to the state 
p2 that results from the choice of q2, for ps is the feed state to 
the final stage. Thus for any choice of g2 the greatest value of 
(6) will be 


i Vip2) — Vips) — Clq2)} + fil pe) 


and vy varying our choice of q2 we shall find the maximum of 
this namely, 


fps) = Max {| V(p2) — Vip C(q2) + filps)}... (7) 


Again this is no great burden since only a few design variables 
have to be simultaneously chosen. Proceeding i in this w ay we 
make the two stages, whose optimal policy and profit we now 
know, the last two of three stages, and in seeking the optimal 
three stage policy vary only the conditions in the first stage 
and always use the optimal two stage policy with respect to 
the state this produces. Thus in ~— 


fr Prsi) = Max | V( pr) —_ V( prea) (qr) )+ fr-i(Pr)} ste 


and the maximization is attained by correct choice of gq. 


To see how vast an economy this can introduce consider a 
six stage process with two design variables, qg, at each stage 
and two variables, p, specify ing the state. Suppose that a search 
over ten values of each of the variables in q is sufficient to deter- 
mine the maximum then if the conditions in all six stages were 
varied simultaneously in the search for the maximum 10!2 
calculations would be needed. By dynamic programming only 
10? calculations are needed in determining the maximum at any 
one stage but f/g—; has to be done for a sufficient range of values 
of the feed state Pe to allow of its being used for the next step 
in Equation (8). Suppose this also requires ten values of cach 
variable in p, then 10* X 10% = 10% calculations will be needed 
at each of the first five stages and so a total of 50,100 calcula- 
tions are needed for all six stages. The undirected search would 
thus take at least ten million times longer than the dynamic 
program. ‘This is the ratio of some 20 years to one minute! 
Of course the unorganized simultaneous variation of all 12 
design variables is about the least intelligent way of tackling 
the problem and it will immediately be suggested that a steepest 
ascent would be much better. 
it reduced the number of calculations to that of the dynamic 


program it would still produce far less information. — For 
the dynamic program produces the optimal r-stage, policies 
r = 1, 2,...R—1 for a whole range of feed states, whereas 


a direct search can at best only provide them for a few, rather 
arbitrary states, and in an intelligent design all such information 
is needed. 

To conclude this general introduction let us see how the 
results of a dynamic programming calculation can be presented. 
Table 1 is constructed with K major sections and, given the 
feed state Peri, We enter the table in section K at this value. 
Corresponding to it we can read of from the succeeding columns 
the optimal policy for stage R, the resulting state and ‘the maxi- 
mum value of the profit function. But the remaining stages use 
the optimal (K—1)-stage policy with respect to their feed 
state pz, so that the conditions in stage K—1 can be found by 
entering section K—1| at pg as shown. This course is pursued 
all the. way up the table and the complete optimal policy 
extracted, Notice that it could have been found equally casily 
for any fewer number of stages. We also notice that the 
algorithm (8) 1s ideally suited to the digital computer for it is 
completely recursive and makes minimal demands on storage 
capacity. “‘0’’ denotes the relevant entry in the table 


The Adiabatic Sequence of Stirred Tanks 


It 1s well known that a revel ible « xothermic reaction goes 
most rapidly if the temperature is made to decrease as the 


reaction proceeds. At any composition there is an optimal 
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[his is indeed true but even if 


























TABLE 1 
No. Feed Optimal | Product | Optimal profit 
of stages state policy state function 
1 | po qi pi fi 
eee 
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os ———0) 
R-1 PR gr-1~ Pr-i fr-1 
eset eeecaen ative 
R Pr+i ‘oR~ sa Pr fr 
“oe ~ 
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temperature at which the reaction rate is greatest and this is 
the temperature that should be used%:®. A rather simple 
method of attaining this by bypassing cold feed suggests itself 
and is illustrated in Figure 2. Of the total flow q a fraction 
Arq 1s preheated to a temperature Tg; and fed to reactor R. 
The product is now mixed with a bypass flow of the unreacted 
feed at its original temperature T, to make up a flow of Ag—~ig 
to reactor (R—1). This mixing with cold feed takes place 
before (or, since they are well mixed, equivalently in) each 
reactor, until finally all of the original flow passes to the last 
reactor (reactor 1) and emerges as product. Each stage is 
insulated and the only heat added to the system is in the pre- 
heater so that it should be possible to achieve an optimal sequence 
of temperatures by rightly choosing Tesi, Arg, Ag-1,..-A% 
subject to Ag S Api S € Av F< A, = 1. Bur these 
quantities will have to be very nicely chosen for too little 
bypass will not decrease the temperature sufficiently, whereas 
too much will retard the reaction. 
n 
If the single reaction }) a:4; = 0 is taking place between 
1 
the » chemical species A), 
any time can be written 


A,, the concentration c; of A; at 


Ci = Cio + QiC.... ee) 


Here c is the extent of reaction and ¢;, the concentration of A; 
in the original feed. ‘The state of the process stream in and 
leaving reactor r can thus be described by two quantities, the 
extent of reaction, ¢,, and the temperature, T,. The cold 
unreacted feed is in the state c = 0, 7 = 7). The problem 
we will consider is that of obtaining the greatest extent of 
reaction in the final product (i.¢, maximizing ¢,) in a sequence of 


tanks of equal volume, V. Now since A, and Cr+, = O we 
can write this 
R 
a = Me Arti CRH = > (Nvle = Avis 6 . (10) 
1 


where Agi Ax as there is no bypassing before the first stage. 


(CoM) 





Figure 2—The sequence of adiabatic stirred tanks with 
bypassing. 
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The rate of reaction, namely the rate of change of ¢ by 
reaction, is a known function of the concentrations and temper- 
ature and so, by (9), of ¢ and T; we write it as r(c,T). A mass 
balance of A; over reactor r and a heat balance yield the equations 


Negeri — Nvt1gCirst — (Ar — Avtt)GCio = iVr(c,, T;).... (11) 
and 

NT, — NenigTr41 — (Ae — Apast)\GT. = HVr(e,, T,),.... (12) 
whiere 


H = heat of reaction/mean molar heat of reagents 


By means of Equation (9) all the Equations (11) i = 1,...” 
become a single equation 
AC = Art iCr+1 - Vr(c,, T,)/q sae (13) 


If we make a slight ~~ of temperature variable and write 

= (T —T,)/H.. (14) 
and Res) = (CV /Oe.T). os cess: (15) 
then (12) and (13) can be written 


Vt, oi Ar+ilrt1 = AC as Ar41Cr+1 = R(¢,, Be Je . (16) 


It is necessary now to distinguish between the first stage 
reactor R with its associated preheater and the remaining stages 
for the first is not truly adiabatic and involves a cost of heating 
that the others do not. Let 

r 
SA Cotter) = Max >> (A,¢, — Ne410s41), 7 ={....R—I, -€27) 
1 
be the maximum value of the profit function attained in the r 
adiabatic stages terminating with reactor 1. For this maximiza- 
tion the bypass ratios A,+1:, A,,. . . Az have to be chosen optimally. 
Applying the principle of ic embodied in Equation (8) 


Fil Cr+ tytras) = Max 1 X,c, = Ars +1 + f yey 4 2) | -. (18) 


where now only X,+: has to be chosen since \, is known from 
the (r—1)-stage policy. The choice of \,+: is restricted only by 


OS Ant & Reside oven a (19) 


For the whole process we should take account of the cost of 
preheating to a temperature Tp+i(= Ty + Htg+i/ep). Let this 
be expressed as a fraction of the value of a unit extent of reaction 
by the function /(t,4:). Then the final choice of preheat temper- 
ature will be the value of ty; that maximizes 


ArCR + fr-1 (Cr, tr) — A(trsi).... (20) 


where Arce = Ar (tr — try) = Ri cr, Phe vi aie c36 a8 0 (21) 


Equations (16), (18) and (19) can be solved for r = 1,2,...R 
to give Xa, Ax and then Equations (20) and (21) can be solved 
for tei:. This is a simple enough matter on the computer and 
little more could be said were it not that a rather simple graphical 
construction illuminates the whole problem. 
We observe first that by Equation (16) we can write 
R 
ec = >} RG, t,), so that we have to maximize the sum of the 
i 


reaction rates at each stage. Now in a plane with coordinates 
cand ¢ the curves of constant R are as shown by the full lines 
in Figure 3, the upper curve corresponding to a smaller value 
of R. Consider first the determination of 


Si(C2,te) = Max(c; = Ax») = Max R(e,,t).. ee (22) 


The point (4, f1) = (Axt2, Aef2) lies on the line from the origin 
through (¢2, f2) and represents the state of the stream after 
mixing with a fraction (1 d») of cold unreacted feed. Indeed 
by Equation (16) which can be written 


ti a= ey R(ei,t) (23) 


we see that it is the feed state to reactor |. We also see by this 
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Figure 3—Optimal conditions for a single stage. 


equation that in the ¢, t-plane the feed state (@, f) lies below 
and to the left of the product state (¢, 1) by a distance R(¢, th). 
We may thus draw a set of loci of &, # for which the reaction 
rate R(c¢, t) 1s constant; these are shewn as broken lines in 
Figure 3. Now for any given (¢2, tz) it is evident that (4, %) 
should be chosen to be the point where one of the broken 
curves is tangent to the ray through (¢2, t2), for this will give 
the largest value of R and so by (22) of the increased yield. 
Evidently this point (é;, 7) is the same for all points (¢2, tz) on 
the same ray (i.¢e., having the same value of the ratio ¢2/t2), 
though the value of A» ae to attain it varies with the 
point (¢2, fs), for A» ¢>. However it is easily possible to 
draw the locus, T,, of F the tangent points (@, 4) and the locus, 
[;, of corresponding points (¢, t) given by Equation (23). 
With these two loci the construction is immediate. Given any 
state (C2, fe) we join this point to the origin and draw a line of 
unit slope from the intersection of this ray with P, tol. Asis 
the ratio of the distance from the origin to the intersection with 
TI’, to the distance from the origin to (¢2, tz) and f, is the difference 
in ordinate between T, and I). Evidently fi(e, t2) = 0 if 
(¢2, te) lies to the left of T,, whilst if (co, t2) is to the right of 
Ts, files, te) = (¢2/te) is a function only of the ratio ¢2/te. 


Consider now two stages for which by (16) and (18) 





fo(c3,ts) = Max} R(c2,t2) + files,te)} 


In looking for this maximum we might as well look for it in 
the easiest way possible, namely along a ray from the origin 
since f; is constant on this. But R is greatest if (ca, t2) is a point 
of tangency of the ray with a curve of constant R. The locus 
I. of optimal states (¢2, fz) can thus be drawn immediately and 
the locus I’, of states representing the optimal feed is obtained 
by taking a point a distance R(¢2, t2)/ Az below and to the left 
ot (¢2, te) tor 


C2 (A; Az)es = ty (As Ne )ts = Ri 2,t2) Ae 


and Xz is known from the one stage policy. Thus for any point 
(¢3, ts) the point (#2, tf) = (Ass ‘Ne, Asls/ Az) lies on Fy and we 
see again how Xg3 and /, can be found graphic: ally. It ts again 
clear that f2(¢s, fs) is zero to the left of fT. and a function only 
of the ratio ¢3/ts to the right. It follows th at in seeking I’; the 
locus of optimal conditions (¢s, ts) we shall be led to the same 
curve as I's, the locus of points of contact of rays with curves 
of the family = constant. 

In this way a system of curves can be drawn Ty, T2,. . . Te-1, 
r,, f.,... Pei, where by coincidence M2, he " are coin- 
cident. Itis illustrated in Figure + tor RK = is the point 
(ca, t4) and this joined to the origin 0 intersecting Is at B. 
BC is the line of unit slope with C on Ts = P and D the 
intersection of DC and [;. We proceed in this way up to G 
on T'\, and have imanedicnchy that f3(c4, ts) is the difference 
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Figure 4—Graphical construction for several stages. 


in ordinates of G and B. The bypass ratios are given by 
de = OF/OE, X3/X2 = OD/OC and d3/A3 = OB/OA. 

For the final stage, reactor R, we will take a simple form 
for the preheating cost and assume that it is proportional to 
the amount of heat supplied, h = aXgte+i. Then by Equation 
(20) we have to maximize 


ArCR + fr (crite) — aArtrs1 


This is a simple direct calculation in which any point (0, te+1) 
is taken and the point (cg, fg) on the line of unit slope through 
it which satisfies Equation (21) is found by trial and error. 
Then the quantity Cr + fr-i/AR can be plotted as a function 
of te+:. For any given a the optimal tr+1 is the point where this 
curve has slope a. Thus the t-axis could be graduated in values 
of a and the curve I'g drawn and the construction would be 
immediate. 

Needless to say the whole calculation would be a matter 
of seconds on a digital computer, which would be very suitable 
for controlling the reactor. Suppose that the preheater is a heat 
exchanger dependent on some other part of the factory and so 
liable to uncontrolled variation. Then a measurement of te+1 
could be fed to a computer and the computer activate the valve 
adjustment mechanisms to keep the bypass rates optimal. It 
would not be much harder, though the calculation would have 
to be rather differently formulated, to control for variations in 
the total flow rate or in the basic feed temperature. If the 
function r(c, 7) were sufficiently well established it might be 
that the computer should be used to calculate a set of tables 
from which these adjustments would be made by hand. 


The General Bypass Design Problem 


The problem we have used in the previous paragraph to 
exhibit the method of dynamic programming is a particularly 
simple one. We now give a general algorithm which allows us 
to determine the size and rate of heat removal, as well as the 
bypass rate, at each state which will maximize a very general 
expression for the profit. 


n 

Let m independent simultaneous reactions }> aj4; = 0, 
i=1 

j = 1,...m, be taking place between the n species, Aj... An. 


Then in place of Equation (9) we have to write for the con- 
centration of A, 


Go = Gio + ze Gor, 8 @ 1. uBio cs. . .(24) 
where c’ is the extent of the 7” reaction and ¢,, the concentration 


of A; in the feed. The time rate of change of ¢ by reaction 
is 7, a known function of composition and temperature and so 
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by (24) a known function of the extents c',...c™ and T. Then 


a mass balance for A; gives 
™m 
Ar+1)9Cio = V, pa ayri(c!,.. .c*,T) 
j=1 


AQCir sa Av419Cinr+1 — (A, 


or by substitution from (24) 
q be ais {Are1 Oar — Ayci, + O,ri(c,,...0%,,T,)} = 0,..... (25) 
=1 


where 6, = V,/q and c’, is the extent of reaction in reactor 1. 
But in saying that the reactions Za;4; = 0 were independent 
we meant that the only set of numbers (,...8, such that 
m 


> a,;8; = 0 is the trivial set Bj = Bo = ... = Ba = O 
i 


It follows therefore that the bracketed expressions in (25) must 
be zero and we have m equations 


Nice Os: — ONG, OTe) Oh iivciiccscccn OOD 
gad... ie 


A heat balance will give an equation 


m 
AT. +9, >> Hr -Q, =0..... (27) 
j=1 


Ari r41 — ArT, — (A, 


where the H; are the heats of reaction divided by the mean 
molar heat of the reactants and Q, is proportional to the rate 
of heat removal from reactor r. 


In the optimal design problem we now have 3K quantities 
to choose for R reactors, 0. . .Az, Qi,...Qpr, A2,... Ag and T+. 
Because of the non-linear way in which T enters the equations 
it is probably more convenient to choose 7),...T, and use 
Equation (27) to calculate Qi,...Qzg. The profit function may 
be constructed quite generally. Let v(c!,...c™) be the value of 
a unit volume of the process stream when the reactions have 
proceeded to the extents ¢',. . .c". Then the gain in value by the 
reaction 1s, 
rts .c™) — 0(0,...0)} = 

q : {r,D(cl,.. 0%) — Negi DWChar,...C™4i)} oe... (28) 


where O(c',...c*) = oc!,...c*) — of0,...0)...........(29) 
The cost of —™, and operating a reactor of volume 
V, = q 6, at temperature 7 * will be a function C(V, rl, Mg) = Cr 
and the cost of preheating to a temperature 7; can be absorbed 
into the expression for this cost in reactor R. We therefore 
have a net profit 


R R 
X pb = u lo [A — bs ®n} — G),.......5. (30) 


where @, denotes O(c',,. . .c™,). 

When the design variables are chosen and the feed state 
Clesty.-.C™e+is Tet is specified Equations (26) and (27) can 
be solved and the net profit Pz calculated. Its maximum value 
is thus a function only of the feed state, 


Sa Onray Dine) = MOk Popiiicvicciciceceni es (31) 


Apply: ing the principle of optimality we are led to the sequence 
of equations 

flCir4t,Tr41) = Max{p, + f-ule,,T,)}, r = 1,...R —1,.....(32) 
where the maximization is by suitable choice of V,, T, and ),+1. 
This maximization will reveal whether or not a bypass is needed 
(A, = A,+1 for no bypass) or if the reactor should be adiabatic 


(Q, = 0). The design variables will be subject to certain 
restrictions probably of the form 
LS Ge oar OR Gy SS ee a  -e \eee | ) 


but this makes the search easier for it limits the region within 
which the maximum must be found. 
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In the final step with feed composition ¢ = 0 T = T, we 
have 
Fr(0,T.) = Max { Pr + Fr-(Cr, Tr)} Sige aire enters (34) 


where pz involves the cost of the arrangements for preheating, 
if these are economic. 


A Limiting Case 


It is interesting to consider the batch reaction as the limiting 
case of an infinite number of infinitesimal stirred tanks. We 
cannot take over the model of section 2 with tanks of equal 
volume, for the actual holding time in stage r is V/A,g. Thus 
letting V0 would lead to infinitesimal time intervals of varying 
length and the limit would not describe the batch reaction. 
Rather we must let the volume V, of the stage vary so that 
the holding time is constant, i.e., V, A,V, where V is the 
volume of : stage 1. Then Equation (16) ‘eonme 


At, = Ap+1lp 41 = A, Cr — NeeiOres = r, Vr(c,,t,)/q Waka ane ak earn (35) 


Defining f,(¢,+1, t-+1) in the same way as before we have the 
same general Equation (18) and the same one stage policy. 
However the equation for two stages now reads 

folcs,ts) = Max{Asce — Ascs + filcs,te)} 
Max{ Az Vr(c2,t2)/q¢ + files,te)} 


If we look for this maximum along a ray of constant ¢2/t2, f; is 
constant as before and Xz is inv ersely proportional to ¢2. This 
time therefore we do not look for the maximum value of r (the 
point of tangency with curves r = constant) but for the maximum 
value of r/c2. But on the ray ¢2 = Pte, r/cz = r(C2, ¢2/B)/ce 
and this has a maximum where 


22 (rcnex/B)\ _ Or | Or 
- z( “Tae ' . 


This defines a fixed locus in the (c, t) plane, which is the pro- 
jection on to the (¢, ¢) plane of the line along which a straight 
line from the origin would touch the three dimensional surface 
r(c, t). The locus of corresponding inlet points 2 = (A3/X2)¢3 
and tz = (A3/A2)ts is now given by 


~~ = Ce => te —_ bo => Vr(€2,t2)/q, 


and since V/q is constant this also is a fixed locus in the ¢, t 
plane. These two loci can be labelled I’; and r, respectiv ely in 
accordance with the preceding example, but this time we see 
that as the construction proceeds all the I’, curves will coincide 
with I’, and all the f, with f,. What is more as V0 these 
two curves will approach one another, so that for a very large 
number of very small stages the system would be described by 
a line zig-zagging between two adjacent curves , and T2. In 
the limit V = 0 these coincide in a single curve I’ whose 
equation is cr, + tr, — r = 0. 

Now the limiting form of E quation (35) can be obtained by 
letting 7 = rV/q be the current time measured from the end of 
the process and 6 = RVq be the total time. Now let V0 and 
R-— so that 6 is constant, then A,, ¢, and ¢, are replaced by 
(7), c(r) and t(r) and 


d d 


Since in a batch reaction the quantity that is directly controllable 
is the rate of addition of cold reagents, namely d), ‘dr, we might 
write the equations 


dX 
— = —4fp,..... : . als: ace era 
dt 
dc Mc 

= — rc. Bi 3 
dr AY t a (37) 
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— = — rc) +—..... , . (38) 
i, 


X is the fraction of the final volume that is there at time 7 and 
since this is increasing as T decreases to zero, u Is positive. 

It would appear then that if the batch process controlled by 
the addition of cold reagent is indeed the limit of the stagewise 
process then it should be induced to proceed along the curve I’. 
Two questions then arise: in the first place how is pu to be 
controlled to keep on I’ and secondly if the state of the system 
is not on I how can it most rapidly be brought there? 


The first question may be quickly answered. If t = y(c) 
is the equation of I’ elimination of r between (37) and (38) and 
integration gives 

{y(c) —c} =A {y(a) —a}..... oa eee 
where ¢ is the value of ¢ when A = 


= Wievy(c)] {U — v'(o/Iy(c) — ev'(o]}.. (ere ee 


The required value of u is 


The second question requires rather a different approach. 
Let c(@), t(@) and X(@) be given and let us ask for the control 
action u(r), 8 > r > O, that will maximize the net conversion 
\(0)c(0) — A(@)c(O). Let this maximum conversion be 


8 
f(u,v,w,0) = Max | ANGI ERE oe occ 
where ¢, t and X satisfy the Equations (36)-(38) and the initial 
conditions 


c(@) = uw, 40) = vo, XO) = w...... 62.25. (4D) 


Now if ¢ is a positive quantity less than @ 


6 6—a 
f(u,v,w,0) = Max ! | Ar dr + | Ar ir] 
O-c a 


The maximization is still by choice of 4, @ 2 r+ 2 O, but if in 
the second integral we do not make an optimal choice of yu, 
6 — oa 2 Tr 2 O, with respect to the initial conditions ¢(@ — a), 
t(@ — a), \(@ — a), we shall certainly not attain the maximum. 
Thus Bellman’s principle of optimality leads us to write 


f(u,v,w,0) = 
@ 
Max | | Ardr + f[c(@—a), (A—c), Noo), 0-2 | ~~ (SS) 


b-o 


and in this maximization yu has only to be chosen in the interval 
6 > 7 > 6 —a. Inthe limit ¢e—0 only u(4) has to be chosen 
and the limit of Equation (43) is the partial differential equation 


fy = Max(r(u,v) {w+ fu + fe} + ufo — (Ufa + rfe)/wh].....(44) 


Here the maximization is by choice of u = 
is obvious. If 


u(8) and this choice 


A = fo — (tha + tf,)/w..... + t's} 


is positive w should be as large as possible; if A is negative u 
should be as small as possible. There will be a limit to the rate 
at which cold reagents can be added so that 


OS pS Fee. ....(46) 


vif A > 0. 


This type of control in which the controlled variable is at 
one or the other of its extreme positions has been variously 
sty led on-off, relay or bang-bang. A is evidently a discriminating 
function saying when adding more reagent is worthw hile 
(A > 0) and when it is not (A <0). The nature of the situation 
is such that when it is worthwhile to add more it should be 
added at the greatest possible rate and that when it is not worth- 
while it should be cut off altogether. The curve along which 
A = 0 is a switching boundary between these two extremes 
and on it «can take any value. We should be inclined to identify 


then up = Oif A < Oand yw = 
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0 D t 


Figure 5—Optimal control of batch reaction by the addi- 
tion of cold reagents. 


this switching boundary with the curve I’ obtained as the limit 
of the stagewise process. This identification can be established 
from Equation (4+) by considering the characteristic equations 
of this partial differential equation. We find first of all that the 
characteristics are reaction paths, so that I’ is a characteristic. 
Then by differentiating the discriminant function A along a 
characteristic we find 


dA/d@ = — {1 + (fu + fo)/w} (ur, + tr, — 1), 


so that A is constant on I’ since its equation is ur, + Ur = 7. 
But this constant can only be zero for were it not « would be 
either zero or v and the reaction path would not stay on I’. 

We should also consider the possibility that the value of 
given by Equation (40) exceeds the upper bound vy on [. In 
Figure 5 the curve ABC is T and on AB the value of given by 
(40) is less than y but on BC it is greater. Let BD be the reaction 
path through B on which » = », then in this case the switching 
boundary should be ABD. The optimal policy can be stated as 
follows: if the state of the system is represented by a point 
lying in the sector OABD the reaction should be allowed with- 
out addition of cold reagents: if the state of the system lies to 
the right of ABD or on the boundary BD cold reagents should 
be added as rapidly as possible: for states on AB the rate of 
addition is given by Equation (40). 
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i 
J 
R 
y 


Nomenclature 


A; = chemical species, 1 = 1,...n 

C(q) = cost of operating under conditions q 

r = extent of a single reaction 

G = concentration of A; 

Cio = initial concentration of A; 

d = extent of j reaction 

fr(pr+i) = maximum profit from R stages with feed state pry; 
H = heat of reaction/mean molar heat of reagents 
h(t) = cost of preheating to temperature ¢ 

Ee = profit function 

b, = state vector for stream leaving stage 

Q, = rate of heat removal from stage r 

= control variable vector for stage r 

qg = flow rate 

R = number of reactors or stages 

R(c,t) = (V/q)r(c,t) 

r(c,t) = rate of reaction in the case of a single reaction 
ri = rate of j reaction. 

S(pr41; dr) = transformation of state effected in stage r 

T = temperature 

is = temperature of cold reagents 

ig.” = lower and upper bounds on temperature 

t = reduced temperature, c,(T — T,)/(— AH) 
u = initial value of c 

V, = volume of reactor r 

v* = upper bound on reactor volume 

V(p) = value of stream in state p 

v = initial value of ¢ 

v(c?) = value of stream when reactions are at extents c/ 
w = initial value of » 

Qj, Qi; = stoichiometric coefficients 

| = proportionality constant 

¥(c) = equation of the curve 

A = uf, + of, — wh» 

0 = total time of process 

6, = holding time of r reactor 


ns = amount of reagent present 

r, = fraction of total flow through stage r 
bu = rate of addition of reagents 

v = upper bound of pu 

T = time from the end of the process 

Ss 


uffix Notation: 


value belonging to 7 chemical species 

value associated with 7 reaction 
value in stage R, the first stage of an R-stage process 
value in stage r, the r stage from the end of the 
process 
value in stage 1, the last stage of the process 


ll 


1 
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The Effect of Intraparticle Temperature 
Distribution on the Catalytic 
Effectiveness Factor of a Porous Catalyst’ 


T. AKEHATA*, S. NAMKOONG?, H. KUBOTA? and 
M. SHINDO? 


The effect of intraparticle temperature distribu- 
tion on the catalytic effectiveness factor is derived, 
and its magnitude is estimated using an approximate 
solution. These calculations show that for the several 
cases examined the term containing the effect of 
temperature is less than 10% of that due to the 
concentration effect. 


I reaction kinetic studies dealing with porous catalysts, the 
effect of intraparticle diffusion on reaction rate is taken into 
account by “‘a catalytic effectiveness factor” which was first 
introduced by Thiele“. Methods of calculating this factor have 
been presented by Thiele” and Wagner“ for first order or 
pseudo first order reactions, and by Kubota and Shindo et al” 
for other than first order reactions. In all these cases, however, 
the temperature distribution within the catalyst pellet has been 
considered uniform. 

Recently Prater‘ has pointed out that the temperature and 
concentration within a pellet are simply related, and that occa- 
sionally the temperature difference between the centre and the 
outer catalyst surface may be rather large. Methods of calculat- 
ing the effect of intraparticle temperature distribution on the 
catalytic effectiveness factor are presente: d, along with illustrative 

calculations for a few commercially important reactions. 


Fundamental Equations and their Solutions 


Assuming a steady state, and that a catalyst having any 
shape can be replaced by an equivalent spherical form, the 
concentration distribution of a particular component within 
the pellet may be obtained from the solution of the following 


| ’ 
D ne ; + v,(¢ t) = 0 (1) 
d) ? d 


and the temperature distribution from the solution of 


Ma +e bene (2) 
A— v,(C, = bebiinis see 
dr? r dr 
with boundary conditions 
dc/dr = dt/dr = 0 atr = 0 
cmc, t=, BCP Og i2 eo ose 63036 (3) 


where v; is reaction rate, defined as moles increase of the 
particular component per unit time per unit volume of catalyst 
pellet. 


1Manuscript received June 30; accepted February 10, 1961. 
2Tokyo Institute of Technology, Tokyo, Japan. 
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Equations (1) and (2) may be combined to relate concentration 
and temperature, as indicated by Prater ®. 


DO 
-(? SV eer +4 = g(c)... ’ (4)* 
Vi 
Substitution of Equation (4) into Equation (1) yields 
d% 2 dc 
o{ ss + =) + ale, ey 20 Oe i. so 
dr? r dr 


Numerical Solution: [Equation (5) may be solved by a numerical 
method as indicated by Kubota and ‘Shindo® even when the 
reaction rate is not first order, or is given in graphical or in 
tabular form. The catalytic effectiveness factor defined as 


(4/2 le, ele) d 
4trrv;|c, g(c)| dr 
Jo 


iba = 6 
, (1/6)md,3 vic, (6) 


is found as a function of ¢, and 4,. 


Approximate Solution: The numerical method above is rather 
tedious, so that an approximate method simplifying the calcula- 
tions is desirable. 

Consider first the case when the temperature within the 
pellet is uniform. Equation (1) may be converted into Equation 
(7), by assuming a Taylor expansion of the reaction rate term 
using only the first order term. Thus, 


d%  2de 
o( =~ ) + vis + 0,(e—c,.) = 0..... satel 
dr? r dr 


vis is the reaction rate at the catalyst surface, and v,’ = dv;/de 
at ¢ = ¢,. This approximation of reaction rate assumes a pseudo 
first order reaction, but differs from the first order approximation 
given by Wagner’ and Hoogschagen®. Schematically, the 
two approximations are illustrated in Figure 1. 

Equation (7) may be solved analytically, as Thiele” has 
done, and the catalytic effectiveness factor is given by 


3 ( 1 =) 
E; = - ——}..... 
m \tanh m m 


*Addition of Equation (1) multiplied by (Q/A;) and Equation (2) 
multiplied by (—1/A;) produces 


d*(rf)/dr? = 0 


where f = (D;Q/A,)c —t. From the boundary conditions, f should 
be equal to (D;Q/A;)c,—t,, and Equation (4) is thus obtained. 
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APPROXIMATION 
BY WAGNER AND HOOGSCHAGEN 


APPROXIMATION 
BY AUTHORS 
ie: 





Figure 1—Approximation of general reaction rate to a first 
order equation. 


where 


~ 


d, Us 
m= — — Givcaieseanuee (9) 
2 D; 
m is regarded as a modified Thiele modulus. 

A comparison of calculated effectiveness factors by the 
numerical method and the approximate method has been made 
for the ammonia synthesis reaction. The result is shown in 
Figure 2, where the relation between the value of Ey’ calculated 
by the approximate method and Fy calculated numerically is 
given by curve B, which is fairly close to the line A expressing 
E, = Ey’. 

Assuming now this pseudo first order approximation, consider 
the effect of temperature distribution within the pellet. 


dc 2 dc 
D; - + + Vis + %5'(C — Cs) + %5/(t — ts) = O.....(10) 
dr? r dr 


TABLE 


CALCULATED EFFECTIVENESS FACTORS 
Assumption; d, = 











1.0 oe | ee 
100 atm. 300atm. 


l 450°C oO e 
500°C @ * 
550°C e@ 


0.8 














0.6 
Er 
0.4 
0.2 
cae 
0 0.2 0.4 05 0.8 1.0 


Figure 2—Relation between E,, calculated numerically and 
E,’, caleulated by the approximate method for ammonia 
synthesis reaction. 


where v;, is the reaction rate at the outer catalyst surface, 


Us’ = (dv,;/dc). = ¢;, 
Ves’ = (dv;/dt), i. 


Substitution of Equation (4) into equation (10) leads to Equation 
(11). 


d’¢ ie 2 dc 1 Vis ‘i v:," 4 %’Q ( 0 
un Sa D, i, CMD icc csasts (11) 


Equation (11) is a first order linear equation with respect to 
c and easily solved. The catalytic effectiveness factor can be 


calculated from Equation (8), if # is modified to 


aa -(% +) Diskette re 
1 


FOR A FEw INDUSTRIAL REACTIONS 
6mm, € = 0.5 




















| A | | 
Reaction Reaction Q D; | nal 8 | Ov,’ | 
Svstem Conditions (cal./mol.) | (cm.?/sec.) | - = 4 * e 
vste oO oO cal./mo (cm.?/sec eye D, i, | m orm Ey 
| i i a wo near nna on snipe cae = a 
ammonia synthesis inlet Hy. 75%, Ne. | 12,800 | 1.02 X 10-8 | 4.07 xX 10-3 | —39,200 | 242 59.2 0.054 
on industrial iron 25% 500°C., 300 atm 
catalyst® Zvuz s=0 
SO; oxidation on inlet SOs; 8%, 23,200 $.53 < 10-* | 762 < 10-* —81.7 6.95 26 | Oi 
industrial vanadium 12%, No; 80% 
catalyst | 450°C., 1 atm 
xSO» 5=0 | | 
| 
naphthalene | inlet CioHs; 0.9%, 520,000** > x 10-* 5 X 10-* —6,080 | 435 | 22.5 0.13 
oxidation on air; 99.1% 350°C., 
vanadium catalyst™ 1 atm xCjHs 5=0 | 
; 
ethylene oxidation inlet CoH 4; 3 %, air; 32.9 X 108N > x 10°? | —62 2.5 | 23 | O65 
on silver catalyst@” = 0% ae atm +3 337. 28 
SH, = 0.5 | x«(1—N)*** | 
CO conversion” | inlet CO; 11.7% 10,000 5.5 X 10-? 8 X 10-4 —7.9 0.28 0.83 | 0,92 
| 350°C., 1 atm | | 
xCO;=0 | | | 
| | | 
Cumene cracking on | 510°C., 1 atm | — 23,700 | 2.0 X 10-3 8 X 10-4 —1,470 —52.6 | 32.7 | 0.092 
silica-alumina xCs= 
a lll | 
assumed ailieos : a a ‘ 
Nomenclature; Z; mole fraction, x; conversion, 
** No side reaction is assumed, i.e., Cyp$Hs + 4.5 hy Pg CsH,03; + 2H.0 + 2C0, + 520,000 cal. 
a : C.H, + 0.5 02 = C2H,0 + 32,900 cal. 
***T yy (eR a cine s ‘eed as: 4 2 ori 4 
Iwo parallel reactions proceed a {eH +3 0. = 2CO. + 2H.0'+ 337,280 cal. 


and selectivity N = k,/(ki + ke) 
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If the experimental reaction rate value v;* includes the effect 
of mass and heat transfer within the pellet, the effectiveness 
factor can be obtained from the relationship + between E; and m* 
which is given by Wagner and Hoogschagen‘®, where m* is 


detined as 
De oe sain (Te eee ; 
m 2 (% + >, ea chee 


In order to calculate the effectiveness factor, values of ; 
and D; must be known. Few reliable data are available, however, 
so that for these illustrative calculations \; and D; are assumed 
from the works of Hoogschagen® and Sehr, and others. 
Calculated results for a few reactions at particular conditions 
are summarized in Table 1. It may be noted that the term 
containing the temperature effect Qv,’/A; is less than 10% of 
the concentration term v,;’/D;. Thus it is concluded that the 
effect of temperature distribution within a catalyst pellet is not 
important. 


SeaeecORe 


Illustration 


Us* = Ey %s’, Ves’ = Ey Ms’, 


then 


dy Qs’ Ves’ 
a = 
m 5 V2 i, + D, ) 


m*? = Ey m,? 


Therefore 


m is expressed by a function of Ey from Equation (8), so that a 
relation between m* and E; is obtained. 


+As v;* = Ey V;, 


Nomenclature 

¢ = Concentration of a particular component 
D; = Effective diffusion coefficient in catalyst 
d, = Effective particle diameter 

E,,E;’ = Catalytic effectiveness factor 

m = Modified Thiele modulus 

Q = Heat of reaction 

r = Distance from centre of spherical catalyst 
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t = Temperature 
t = Temperature of surface of catalyst particle 
v; = Reaction rate, moles increase per unit time per unit 
volume of catalyst particle ( = v/(1 — €), where v is 
reaction rate for packed bed of catalyst, and € is void 
fraction of catalyst bed, 
, dv; , 7 ss 
Vs = 7c neglecting the effect of temperature distribution 
dc 
5 _ (00; y 
=~ ee 
_ fax\ | ; 
- — 
Vj = Thermal conductivity of catalyst particle 
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A Design Parameter for 
Multicomponent Tray Design Estimates’ 


ALFRED J. SUROWIEC? 


This report develops computational techniques 
that are useful in multicomponent distillation. A 
design parameter, (dn/0N),, is introduced which is 
numerically equal to the number of theoretical plates 
required to do the work of one total reflux tray at 
the bottom of the enriching section. For a minimum 
sized tower, (9n/QN), is less than 3.5. Economic con- 
ditions may permit values as high as 6. 


The ratio of the number of theoretical trays, n, 
to the minimum number, N, for the combined enrich- 
ing and stripping sections, including reboiler and 
condenser, and the ratio of the reflux leaving the 
enriching section, L,, to the minimum reflux, L,,, can 
be found from: 


2.303(0n/ON), log iol (On/d N)s )s 


/N = 
(On/ON), =. 


7 (On/ON), 
piatesdeets are 


A total of 36 sets of minimum reflux data in multi- 
component systems were correlated with an average 
error of —4% in the ratio n/N by these equations. 


Since (@n/dN), is related to the relative volatility 
and the compositions of the counterflowing liquid and 
vapor streams at the bottom of the enriching section, 
the reflux L, can be found from a feed tray balance. 
A procedure is provided for evaluating L, after a 
value for the design parameter has been selected. 


hen a designer has selected a reflux ratio and made a tray 

calculation, he frequently has little idea of the kind of 
selection he has made and may be forced into an involved 
economic balance in order to find out. Even when such a 
balance has been made, he still has little idea where he stands 
in relation to the distillation process itself. A logical measure 
of the desirability of a design choice is to compare it with the 
minimum number of trays and the minimum reflux for the same 
separation. A more convenient and flexible measure is to intro- 
duce a design parameter which gives the number of theoretical 
operating trays required to do the work of one total reflux tray 
at any point in the tower. 


Basic Equations 1 compare the number of operating trays 
with the number of total reflux trays it is useful to place the 
operating line equations in the same form as the equilibrium 
relationship by introducing a composition ratio ¢, such that“ 


| = OnXn ces esha 6 OREO woe we O.® (1) 


iManuscript received October 16, 1960; accepted March 15, 1961. 
2Chemical Engineer, 78 Midland Boulevard, Maplewood, N.J., U.S.A. 
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Then, along with the equilibrium relationship 
Vu = Gude. SE er eT ee 62) 
there results for a single tray 
Val Vots = Qe/@nx:.. 
In logarithmic form 
log Y, — log Yni1 = log(an/@n)..............(4) 


For differential numbers of trays, the following procedure i is 
a short cut to the lengthy derivation of E quation (71) in 
Appendix B. Differential trays are discussed under Mathematical 
Conventions in Appendix C. 


Using Taylor’s series for the finite difference in Equation (4), 





© )Ad*log Y << 
log(a/)x - ys a a cig Meas (5) 
Differentiating, 
d 2 — )tkd*log Vu+t 
lon a/d)n = — x! tae! (6) 


Combining, 
d d : 
log(a/®)» = > 3 log(a/@)n + log Vn pr 
dn dn 


— (—)K A — k/2)d*log Vn+1 


aoe 7) 
<= k! dn* ( 

Rearranging and neglecting third order and higher terms, 

log(a/d)n dn + 3d log(a/@), + d log Vasu. = 0........(8), (71) 


At total reflux, @ is equal to unity and the minimum number of 
trays is given by 


logaydN + 3d logay + dlogYyy,1 = 0.........(9) 


Equations (8) and (9) were found to have an error less than 
4% in the number of trays. (See discussion of Equation (76) 
in Appendix B.) 

The component distribution at total reflux is different from 
that at some other reflux ratio. However, for the keys, ay 
can be taken equal to a, when Vy+1 equals V,+41. Equations (8) 
and (9) then yield 


logan . } d logdn 


~ Jog(a/d)n z * log(a/)n 


Equation (10) holds for any pair of components where there is 
little difference between a, and ay. In what follows it will be 
written only for the key components. 

Equation (10) shows that the differential of 7 has been written 
as a function of N and 4, i.e., 


= (On/ON)dN + (On/db)db............ (11) 





ney eee . (10) 
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Inspection shows that 
logan 


(dn/ON) = josie /6). eT PE Cre (12) 


Rearranging, 
logan ae (On/ON) 


logdn (On/ON) —1 








Differentiating and neglecting a term in d loga which is not 
justified by subsequent approximate methods of integration, 
logan 


d logdn = (On /ON) QGRsGN). <0. 6.0063 0<€44) 


Making the substitutions in Equation (10) 
= (On/ON)dN + 43d log(On/ON)...........(15) 


Integration of Equations (8) and (15) can be accomplished 
by taking logY,,—logY,+: linear in logY,+:°, and by setting 
1/(0n/ON) proportional to N. Then, 


A log( Y,/ Y, nti) 


o log log( Yn/ Yn+1) 








= Alog Yas: + $A log( Vn/ Yuri). (16) 


and, 
i 
. ‘ba ! 
An.- — : - = AN — 3A (am) ers (17) 
. we cssam) 
where 
ea ae 


"A log logay 


At first glance, there would appear to be no particular 
advantage of E \quation (17) over ( (16). However, the differential 
coefficient defined by Equation (12) is a useful design parameter 
which gives the ratio of the operating trays to the total refiux 
trays at any point in the tower. The linear approximation used 
in arriving at Equation (16) was found to yield accuracies of 
from 5 to 10% in the number of trays. The equations are 
applicable to the stripping or enriching sections but require a 
feed plate match before they can be applied to the overall tower. 


Feed Plate Match A feed plate match can be made by re- 
ducing the multicomponent system to an equivalent binary 
system and by introducing feed at the intersection of the 
operating lines. In the new basis, the operating line in the 
enriching section takes the form 
Vats (m + Yh) n+1 (Y/1 + V)sa4i = 

L(t +X) (X/L AX jn + D(x1+xn)v (X/14+X)p; 


RG ga ee aga cc (19) 

where 
(Y/1 + VY)atr = (OX/1 + OX), (20) 
(¥/Al + ¥), = (aX/1 + OX)o........0..+ 428) 


Similar expressions can be written for the stripping section for 


+2 (f—!). 


For a partially vaporized feed, a component balance gives 


F(z: + 2) (Z/1 + Z) = 


Fi(xi-+txn)e (X/14+-X)e + Frintya)e(V/1+V)e.......(22) 

Woe Wa iicccdin vo dukiwcevsxns (23) 

The operating lines for stripping and enriching are coupled by 

Ly—-1 (41 + Xela — Lex + xnde = Fr(xi + Xa)p...... «so (S4e) 

Vers (a + Yess — VA + yy = Fel + ynde-.-.-.-. (24b) 
The intersection of the operating lines is defined by 

ee I ee a Soot (25) 


Rearranging Equations (24a & b) with the aid of a component 
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balance and allowing ¢ or f—1 to approach the intersection, 


—Fixi + xa)e a + X)p — (X/1 + = 





Fim + ye | (Y/1 + Yor — (Y/1 + V)sas 
LAx1+Xh )e (X/1+X), — (X/ 1+X), 
ee EF ere ...(26a) 
Vers(ntynderr | (V/1+V a1 — (V/A+ Ves bal 
Ly-(Xi Xn )y—-1 (X/1+X), -— (X/1 +X yt 
Vi ni +ya dy (Y, 1+V),+1 ~ (VA+Y¥)y ], 








(26b) 


The right hand sides of Equations (26a & b) tend toward 
unity so that 
_ ~ Fi( ut Xa)e 

Fenty 


_(V/A+V)e — (Z/+Z) 2 


(X42) — (Z/1+2) 


(Y¥/1+¥)u1 — “ (Z/A+2) (27) 
(X/1+X), — (Z/14+Z) ) 


Equation (27) also applies when feed is flashed on an equilib- 
rium feed tray, and 


Xy = Xeur = Xp; Vy = Vous = Vor..........(28) 
For an all liquid feed at its bubble point, 
X, = X, = Xp = Xp; Fy = 0; Z = Xp....... .. (29) 
For an all vapor feed at its dew point, 
Vinx = Your = ¥, = Yr; Fp = 0; Z = Yp..... —e 


Using these feed conventions, and setting (0n/N) equal to 
unity at Xp and Xz, Equation (17) takes the form, for the 
overall tower: (See Appendix D.) 

_ (dn/d N), ont Se/ ON), 


cays (dn/ON), — an 


It is assumed in the above that 
Qe-1 =a, = @ 
so that, for a saturated liquid or vapor feed, 


(On/ON). = (On/ON), = (On/ON),-1 
Figure 1 shows the intersection of the operating lines. 

A certain amount of caution should be exercised in the use 
of Equation (31). Refluxing of excessive light and heavy 
components about the feed tray may cause a reverse fraction- 
ation which can result in the operating lines intersecting on the 
opposite side of the equilibrium curve with ¢ >a. This situation 
frequently occurs at minimum reflux. Theoretically, refluxing 
of heavy non-product components in the fractionation section 
can be carried out until all of the overhead product Components 
have been displaced from the reflux. Then, since V = L + D, 
the product components appear in the vapor stream ‘i the ‘same 
ratio as in the overhead product. Inspection of the particular 
case at hand will show if it is advantageous to assume in (31) 
that all non-product components are absent from the enriching 
and stripping sections when the design parameter is determined 
and then to allow for the refluxing effect around the feed tray, 


Effect of Reflux Ratio When the design parameter is con- 
stant, the ratio of the actual number of trays to the minimum 
number of trays for the overall tower is equal to the parameter 
itself. To obtain, for example, a system requiring three times 
the minimum number of trays with (0”/0N) constant at 3, 

is only necessary to draw in the total reflux steps on a McC abe. 
Thiele diagram in the equivalent binary system, and then, 
between the first and second steps on the 45° diagonal, begin 
two more total reflux constructions. The diagram is now 
composed of three interlacing total reflux steps with es 
given by their intersections. (Note that X, Y,+3-) The 
operating line is now a curve drawn nearly one a tind the way 
from the equilibrium curve to the 45° diagonal, as in Figure 2. 
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(Y¥/1+Y) 


(X/1 +X) 


Figure 1—Intersection of the operating lines and extrapo- 
lation of minimum reflux in the equivalent binary system. 


This tv e of construction can be repeated for any rational 
value of n/N equal to (0n/0N). The geometry of the problem 
>: = that the distance between the equilibrium curve and the 

° diagonal is broken up into nearly equal parts totalling 
a ON) in number. The distance between the equilibrium 
curve and the 45° diagonal is given by 


(V/A+V)nar — (X/1 AX) ng & (V/V) — (X/1 4X) 
while the distance from the operating line to the 45° line is in 
turn 

(VY/i+V)na1 — (X/14+X), 


The ratio of the two distances is then 


(dn/ON) _ (Y/Y ss — (X/14X ost (32a) 
(@n/ON) —1], 9 (V/V enn —(X/4X)e 


(On/ON) (Y/Ai+Y), — (X/1+X), (321 
_ = ks 8 IZD) 
(On/ON) —1 (Y/1+V)ex1 — (X/14+X)s 
A little manipulation gives 
(On/ON) Qn + Vn: Quit — 1 (33 
= i eae I35a ) 
nil gust Pou eye 
(On/ON) 1+ O,X, @ — 1 (331 
(Qn/ON)—1|, 1+40.X,  —1 


Equations (33a) and (33b) are at best approximations. They 
yield the same value for the design parameter when lines of 
constant @ and a are parallel. For small values of relative 
volatilities they reduce to Equation (13). They will be useful in 
determining the effect of reflux ratio on the design parameter. 


Define L’y, as the reflux that obtains when the operating 
line in the enriching section intersects the equilibrium curve. 
L’y, is not a true minimum reflux but a pseudo minimum reflux 
that is obtained from a geometrical construction with an equilib- 
rium curve that has already been established and not from 
material and equilibrium balances. 

When L’y, intersects the equilibrium curve at Y = Y,y41, 
a condition analagous to an all vapor feed at its dew point, 
Equation (19) yields 


Ln Xi Xb) _(Y/+Y hn — (X/1 +X nt (34 
FS  _................. J+) 
Lal X'1+X's)un (VY/1+V)y41 — (X/14+2X), 

(On/ON) 
om ee (35) 
(On/ON) — 1 : 
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(¥/1+Y) 





(X/14#X) 


Figure 2—McCabe-Thiele diagram in the equivalent binary 
system when (dn/0N) equals 3. 


When L’x,,, intersects the equilibrium curve at X = X,, a 
condition analagous to an all liquid feed at its bubble point, then 


Lal x1+Xh)n _ 
L' unl X'1+%'s) Mtn = 
(X/1+X)p — (Y/1+V)y4s (Y/1+Y hie —(X/1+X), 6) 
(X/1+X)p —(V/1+Y), (Y/1+Y ia os (X/1 +X), 
_(Xp —~ eX. we + On. aL (dn/ON) | (37) 
(Xp — a,X,)(1 + Xp) (On/ON) — 1 


In the fractionation section, components lighter than the 
light key component soon reach constant com ositions®», and 
I 


gp OO ON 5 ns cnn he, ea ig EN AT TE Wa RIES weer 
= (L,/1 n+l )Xin + (D/V,, TEP eve e ete eevee ewes (39) 
and 
Dx; 
LiXin = pat oinnengis ase .(40) 


Kin ; 
(L,,/ Vu+1) 


During the refluxing of heavy components in the fractionation 
section, the heavy key component goes through a maximum. A 
maximum implies a region of constant composition, however 
small. Assume Equation (38) applies. Designating wu as the 
tray at which the maximum occurs and making use of the fact 
that the overhead product is small for the heavy key component, 
Equation (40) reduces to 


Ki “ 2 


| fp | GO ae, <a 2 (41) 
A similar treatment in the stripping section for the compo- 
nents heavier than the heavy key component yields 


Bx 
ithe sacri Me MD ath 2! (42) 


\ Kin 
(L,/ Vint 1) 


while for the light key component at its point of maximum 
concentration 
Eig asl Vesa MOS = Bin Pew es (43) 


Relative volatilities can be substituted in Equation (40) and in 
Equation (42) when they are written for n equal to u and n 
equal to w. 


Correlation of Minimum Reflux [Except for n = ¢, all 
values of the true minimum reflux, Ly,,, are essentially the same 
as the pseudo minimum reflux. Since heavy components are 
fractionated out as ” grows large, and lighter components reach 
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TABLE 1 


COMPARISON OF RESULTS OF EQUATIONS (31) & (44) 











= Test No. Feed Cond. N 
(A) Ex. 2 | BP. | 8.64 
3 60% Vap. es 
5 D.F. ~ 
6 a ee 
7** | a - 
(B) Ta 8% Vap. 26.85 
b s * 
c | | ‘i 
Ila si 25.65 
a 4 +3 
II 72% Vap. 
IV Bue. 18.4 
V | " 10.44 
VI 50% Vap. 13.3 
VII | BP. 29.45 
(C) II B.P. | 49 
B.P. ad 
III u 11.2 
IV-A *** 9.4 
IV-B or 43.2 
VII B.P. 9.0 
VIII | = = 
IX | ‘ 29'0 
xX | 20.2 
XI 2 11.8 
XII a nie 
| XII is 86 
| XIV Ks 13.0 
XV 68°  Vap. 3.4 
XVI B.P.. 19.7 





* Not reported. 

** Reverse fractionation in the stripping section™, 

*** Cold feed. 

+ Rejected. 

(A) Murdoch, P. G., Holland, C. D., Chem. Eng. Progr., 48, 
(B) Brown, G. _ , Martin, H. Z., Trans. Am. Inst. Chem. 
(C) Gilliland, E. R., Ind. Eng. Chem., 32, 1220 (1940). 


constant compositions, the sum of the combined keys obtained 
by the intersection of the operating line and equilibrium curve 
on the equivalent binary McCabe-Thiele diagram differs slightly 
from the result calculated from an equilibrium balance because 
of changes in relative volatility. The component distributions at 
the operating and minimum refluxes are not the same. The 
essential difference at 7 = e between the pseudo minimum reflux 
and the true minimum reflux is the quantity of heavy material 
that enters the enriching section from the feed tray. The sum 
of the combined keys remains nearly the same. F urther, a feed 
tray analysis will reveal if there is any appreciable difference 
between the percentages of combined key s at the true minimum 
reflux and the operating reflux. As a first approximation, 
Equations (35) and (37) can be written as 

(L/L), = (On/ON)/[(On/ON) — 1]; n =e,s.. (44) 
‘Table 1 shows the results obtained with Equations (31) and (44). 
Design Procedure To arrive at a tower design, select a 
value of (0n/ON), for use in Equation (31). From Equation (13), 
obtain @, and then by means of Equations (22) through (30), find 


As = Lxu/LXis 
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a Calculated Actual 
(L/D) L/Lu n/N n/N n/N-~Calc. 

0.937 £:25 2a) 2.00 1.06 
* 2.00 a 1.00 

™ 2.06 - 1.03 

. 2.08 1.04 
. = 1.99 0.996 
532 | 6.08 Rh 1.09 1.02 
we 3.57 1.19 1.47 1.02 
a 2.50 1.29 1.28 1.01 
4.37 5.50 1.14 1.10 1.04 
‘ 3.04 1.27 1.21 1.05 
‘a 2.06 1.47 1.37 1.07 
3.25 a 1.33 1.24 1.07 
1.83 2.06 1.47 bos 1.07 
2.54 2.89 1.34 1.22 1.10 
0.942 3.18 1.28 1.20 1.07 
7.00 1.36 1.77 1.80 0.983 
win 1.92 1.42 1.41 1.01 
b.35 1.84 1.82 1.01 

0. ‘98 2.04 1.52 1.38 1.10 
0.53 ade 1.38 1.15 1.20+ 
0.76 2.63 +.52 1.26 1.21f 
1.25 1.28 2.00 1.94 1.03 
a 1.68 1.67 1532 1.10 
a 4.11 1.22 1.14 1.07 
6.0 2.0 | 1.35 1.39 0.97 
7 3.16 1.21 Edt 1.00 
7.0 2.86 23 1.22 0.99 
vs 5.00 1.10 125 0.96 
2:3 Pane | 1.43 1.34 1.07 
“ | 69 | ee | 2 1.00 
2.26 1.14 | 2.35 | 2.4 0.98 
i 3.78 1.23 1.16 1.06 
‘ 6.00 1.14 1.09 1.05 
1.85 4.86 j ey 1.12 1.09 
0.9 3.33 1.31 1.19 1.09 
0.78 1.15 2.94 2.32 1.27+ 
: 2.56 1.65 1.27 130+ 
7.3 1.21 2.09 zon 0.99 
1.82 1.59 1.45 1.09 

3.03 1.30 uaa 1.06 


254 (1952). 
35, 679 (1939). 


Since FY, and X, are known, the sum of the combined keys, 


I s(Xis tXas) = V4 iy + —_ D(xip+Xap) 


1+Y¥as+1) 


can be readily found from Equation (19), by eliminating the 
vapor term. 


The sum of a distributed component and the heavy key 
component can be obtained in the same manner using a value of 
the design parameter equal to that used for the light and heavy 
key components. Components lighter than the light } key and 
heavier than the heavy key can be found from Equations (40) 
and (42) and a feed tray balance. Sufficient information is now 
available to arrive at L,. (See Appendix E). 


Relative volatilities can be used as a first trial in Equations 
(40) and (42). If a trial value of L, is used, it can be set equal 
to L, in (41). Using the same value of the design parameter 
and the same value of # at the intersection of the operating lines 
when fractionating between the distributed component and the 
heavy key component is equivalent to saying that /N 1s the 
same as it is for the light and heavy key components. If desired, 
the feed tray balance can be reworked and (dn/AN), can be 

calculated and compared with that originally used. 
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Design Conditions In an economic study with the aid of a 
computer, Albright‘® found that reflux ratios of 1.20 to 1.30 
times the minimum were optimum for his conditions, while 
Happel‘” presents a correlation which prescribes reflux ratios 
that are close to the minimum. Substituting these reflux factors 
in Equations (44) and (31) yields the following: 





(L/La)s n/N (On/ON), 
| 
1.2 | 2.14 6.0 
1.25 | 2.00 5.0 
1.3 1.9 4.33 


These values should be compared with those employed in the 
sample problems. 

In this connection it is interesting to note that when ~ 
relative volatility is sufficiently small, then Equation (44), 
(53), predicts the results of Cohen“: An enriching section ‘wd 
a minimum size, or volume, when the reflux ratio 1s twice the 
minimum at all points and when twice the minimum number of 
trays are employed. Going a a step further, an enriching section 
with a constant reflux ratio has a minimum size when: 





an, 4 ES i See eee 
n/N 1 a75 | 1.73 1.65 1.00 
N. | 8.73 ia i ~. 
ite 4 BA 142 1.50 oo 


The combined enriching and stripping sections is a minimum 
when (07/0N), is less than 3.5. 


Discussion 

The design parameter introduced herein has certain advant- 
ages: 

(1) It permits a rapid estimate of the total number of trays 

from the minimum number of trays. 

(2) It permits an estimate of the required reflux ratio from 

a feed tray balance and eliminates the need for a minimum 
reflux calculation. 

(3) It allows the designer to select a set of design conditions 

which, if not the optimum, are none the less desirable. 

The use of the pseudo minimum reflux is based upon the 
geometry of the McCabe-Thiele diagram®) alone. It supposes 
that the equilibrium curve has been determined by the separation 
under consideration. 

Conventional methods are employed in making a feed tray 
balance“. Location of the feed tray is a free choice and for 
computational convenience it is taken at the intersection of the 
operating lines in the equivalent binary system. The question 
of optimum feed tray location is left to the consideration of the 
particular case at hand. For “normal” systems, the optimum 
location occurs at the intersection of the operating lines in the 
equivalent binary system. 


The use of Equation (44) in the optimization studies in the 
Appendix A is, admittedly, an over simplification. Its use was 
felt justified because of the results obtained in Table 1. The 
data of Table 1 could be represented with an average error of 
—4% in the ratio (n/N), by means of E uations (44) and (31). 
For the systems represented in Table 1, the key components 
have relative volatilities of 3.0 or less. Binary systems are not 
included in the data of Table | 

Referring to Figure 1, it should be noted that crossing the 
equilibrium curve during the refluxing of non-product compo- 
nents, particularly at minimum reflux, forces the operating line 
to change direction within a single step. With an infinite 
number of steps at the crossing during minimum reflux, this 
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results in a discontinuity of slope of the operating line. An 
infinity of steps at or near equilibrium at the point of crossing 
suggests that the operating line becomes tangent to the equilib- 
rium curve at this point. Thus, if the intersection of the operating 
lines in the equivalent binary system is known for the minimum 
reflux, the pinch points in the fractionating and stripping sections 
can be found approximately, simply by drawing tangents to the 
equilibrium curve from the point of intersection. 

Use of the linear approximation in the integration of Equa- 
tions (8) and (15) follows from the fact that the operating line 
for enriching can be put into the form 


(Vint) (Y-1)/(¥+1)ag = 
[L(xr+xn) (X —1)/(X +1), + Dr —xa)d 


In terms of logarithms, 


log Vnsi = Ln 21+ h)n log Xn + 2D(x1—Xh)d + 2¢ 
Vins W+Yn n+ 


where 


ws i V-1 1+2k Lal 1++-%u)e -(33 1+2k 
gee. = ae, art aera ; 
Node Peet Ne +17, Var +¥a ast \X +1 b 


When Fn41 

When the operating lines in the equivalent binary system 
are straight, then the slope of the logY,,+: vs. logX, curve goes 
through a maximum in the enriching section and use of the 
linear approximation in the integration of E quations (8) and (15) 
can result in values of 7 somewhat larger than actual. Usually, 
L(x: +*s)» increases with increasing ” for n < e, and then 
decreases as the heavy components begin to appear. If the 
intersection of the operating lines is determined from the 
maximum flows, then the linear approximation will yield values 
of nm less than actual. At constant total molal downflow, the 
maximum flows for the combined keys are given by: 


< 2 and X, > }, then g becomes small. 


A 


Lil xi +X: \a,max = LA XIAN hy + zu Ve+i Vie+1 , nN £e. dick Hob 86 .(45) 
i> 


IIA 


Ly Xt + Xn yaa + 2», Ly-1 Xig-1; n >f—1 
i 


A straight operating line in the equivalent binary sytem 
permits the use of the coordinate transformation of Stoppel® 
in Equation (9) by replacing Vy41 with 


(X/1+X)n — (X/1+X) 


(X/1+X) — (X/1+X)p 


(V/1+¥ ues — (LAFTL) _ 
(V/AI+Y) — (Y/1+¥)n41 


where (X, ¥) and (X, Y) are the left-hand and right-hand inter- 
sections of the operating line and equilibrium curve. 
Intersections of the operating line and the equilibrium curve 
are obtained by assuming a to be constant outside the region of 
immediate interest. Such regions should not include irregular 
variations of volatility such as that which may occur across a 
condenser. Thus, in the enriching section, for example: 


n Ss 


Y=a,X; Y=a,Jd; 1 


IIA 
IIA 


rar Xt +X) 


n+ (+I )ov1 > 


(V¥/1+Z) — (X/1+X)p _ (V/AL+Y) — (X/1+X)p 
(X/14+X) —(X/14+X)p (X/14+X) —(X/14+X)p 





It should be noted that X can be negative, passing from + © 
to — © when X/1 +X is unity. 
If the relative volatility is constant, then ay can be replaced by 


1 + (G@,—1) (X/1+X) 


1 + (an—1) (X/1+X) 





If ay is not constant, then the transformed volatility is calculated 
in the obvious way from known quantities. With suitable 
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assumptions, a final form for 7, analagous to Equation (18), is 
obiained. 


Nomenclature 

B = Bottoms product, moles 

D ~~ = Overhead product, moles 

e = Theoretical plate number at the bottom of the enriching 
section 

f—1 = Theoretical plate number at the top of the stripping 
section 

F = Total feed, moles 

F, = Liquid portion of the feed, moles 

Fy = Vapor portion of the feed, moles 


k = Index of demarkation 


K,, = Equilibrium constant for component 7 at the temperature 
and pressure of plate 

log = Natural logarithm 

Lye = Minimum reflux leaving the bottom tray of the enriching 
section obtained from a feed tray balance, moles 

Ly,» = Minimum reflux leaving plate ” obtained from an equilib- 
rium balance, moles 

Lu; = Minimum reflux leaving the enriching section at the 
intersection of the operating lines given by Equation 
(44), moles 

L, = Reflux leaving plate n, moles 

L, = Reflux leaving the enriching section at the intersection 
of the operating lines, moles 

L'wn = Pseudo minimum reflux leaving plate m obtained by a 
geometrical construction, moles; n <e 

n = Theoretical plate number counting from the top of the 


distillation tower. A partial condenser with condensate 
to reflux is taken as the first plate, and an equilibrium 
reboiler as the last. See Appendix D 
N= Total reflux trays for the key components computed in 
terms of vapor compositions 
(On/ON), = Differential coefficient defined by Equations (11) 
and (12) and evaluated at the intersection of the 
operating lines in Equation (31) 
Value of n at the intersection of the operating lines 


Ss = 

u = Theoretical plate in the enriching section where the 
heavy key component reaches its maximum concentra- 
tion 

Vi+1 = Vapor stream entering ple ate m, moles. 

V.+1 = Vapor entering the enriching section at the intersection 


of the operating lines, moles 
w = Theoretical plate number in the st: ipping section where 
the light key component is at its maxitnum concentration 
= Mole fraction of component 7 in the liquid portion of 


NiF oa 
the feed 
xi, = Mole fraction of component ¢ in the liquid leaving plate 
X = Mole ratio: moles of light key component per mole of 
heavy key in the liquid 
Xz = Mole ratio in the bottoms product 


Xp = Mole ratio in the overhead product 


X; = Mole ratio in the liquid portion of the feed 

X, = Mole ratio in the liquid leaving plate 

X, = Moleratio in the liquid at the intersection of the operat- 
ing lines 

yir = Mole fraction of component 7 in the vapor portion of 
the feed 

Yint1 = Mole fraction of component 7 in the vapor entering 
plate n 

Y = Mole ratio: moles of light key component per mole of 
heavy key in the vapor 

Yr = Mole ratio in the vapor portion of the feed 

Yau = Mole ratio in the vapor entering plate n 

Y.41 = Mole ratio in the vapor at the intersection of the operat- 
a lines 

2; = Mole fraction of component 7 in the total feed 

Z = Moles of light key component per mole of heavy key 


component in the total feed 


Greek Letters 


aij, = Volatility of component 7 relative to the volatility of 
component j at the temperature and pressure of plate n 
A = Final value minus the initial value of the quantity to 
which it is applied over the region of interest concerned 
dijn = Composition ratio defined by Equation (1) 
Superseripts 


Pseudo minimum reflux 
Maximum flows 


” 


ll 


Subscripts 
h = Heavy key component 
! = Light key component 


Minimum reflux 

In general the subscripts / and h are not used unless it is 
necessary to distinguish between components other than the light 
and heavy keys. 


M 


ll 
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APPENDIX A: OPTIMIZATION STUDIES 


Since the reflux loading in a tower determines to a large 
extent the diameter, the size of an enriching section can be taken 
as proportional to the integral 


e 


P = | (L,/D)dn 


0 


To find a minimum sized plant it is necessary to determine the 
conditions under which the integral P is stationary. 


Variable Reflux Case. Inspection of Figure 2 shows that the 
reflux loading Li(xm + Xi») can be gotten in terms of (dn/0N) 
and Y,41, or with Equations (9) and (15), in terms of » and 
(On/ON). Ina region where Equations (40) and (41) are applic- 
able, the percentage of combined keys has a value which depends 
only upon L, when a does not vary. The functional dependence 
of the reflux can then be indicated by 


(L,/D) = L/D |(On/ON), n]..... ..(46) 


According to the Calculus of Variations“, the integral P will be 
stationary when 


O(L,/D) 

A(L,/D) Op . 
0(0n/ON) dn a ee (47) 

_ ddn/ ON) 
Reference to (46) shows that p = — does not appear 

n 
and Equation (47) reduces to 
O(L,/D) 

=O (48) 


0(0n/ON) 
Equation (48) requires that the design parameter remain con- 
stant, hence, 


dn = (On/ON)dN......... (49) 
Rewriting the integral P 
a N. 
= (0n/ON) (L,/D)dN..... me) 
0 


Differentiating with respect to the design parameter and setting 
the result equal to zero 
dP 
d(dn/ON) 
Nef L, N. A(L,/D) 


dN /ON) - N. 5 
D 7  aan/ON)* es. 


ihc 


Making use of Equation (1) and the fact that (x1'/x4')an = Xn+1 
in Equations (35) and (33a) results in 


Lk Qn+1 — 1 Xie 


— = Sn a bisa dicvace iw s(S2) 
L' vin Pn — 1 Xin ( 
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with X17, 2 Xmai. With a sufficiently small and of little varia- 
tion, Equations (52) and (13) yield 


L. ~ (On/ON) 
Lo we = (On/ON) — 1 


.(53) 


Substituting (53) in (51) results in 


(On/ON) rN. 


(On/ON) a ; “ 
: (an 3 V) Te (L Un D)dN = 0..(54) 


(On/ON) —1 


Equation (54) reduces to 
(On/ON) = 2 
and (L/Lu), = (L/L'u)n = 2; n Ce 
Constant Reflux. Setting Z, equal to LZ, in the integral P, 
P = al... ; wi 433) 


Rewriting Equation (55) for » = e with the aid of (44) and (17), 
differentiating with respect to (0n/ON), and setting the result 
equal to zero, 


, s(logm — m + 1) 
N, = 7m 


a: - (On/ON). ......(56) 
(m/m — 1)(m — 1 — 2 logm) 


The solution of (56) yields the following for a minimum sized 
enriching section with constant reflux: 


N. 


(On/ON), 


0.00 1.00 
1.00 2.82 
1.54 3.00 
2.02 3.10 
8.73 | 3.40 
16.10 3.45 
0 3.50 


For the stripping section, P is equal to (n — s) (L; + Fz), 
while for the combined enriching and stripping sections, P is 
given by nL, + (n — s)F,. The overall tower is a minimum 
when Ly;N/Frz is equal to 
(m — 1)? (1/2m) + (N — N,) (logm — m + 1) 


Rate . im 


(On/ON), 
(m/m — 1)(m — 1 — 2 logm) 


In order for P to be a minimum for the overall tower, (0n/ON), 
has a maximum value of 3.5. 


APPENDIX B: DERIVATION OF EQUATIONS 


The basic differential equations of multicomponent distillation 
have been arrived at by resorting to a standard technique and 
employing Taylor’s formula. To find the differential change in 
n, expand v at the point log Y,4:, logX,, and loga,, by the Taylor 
formula, retaining only first order terms. Then, 
dn = (On/0 logXn)y,ad logXn + (On/0 log Vny1)x,0 dlog Vn41. .(57) 

+ (0n/0 loga,)x.y d loga, 
Making the substitutions 


d logX, = d log Yni1 — d logdn........ (58) 


and 


d loga, = d log(a/@), + d logd, (59) 


Equation (57) becomes 


dn = [(On/0 logXn)y.q + (On/0 log Vnii1) xa] d log Vnsi......(60) 
+ [O0n/0d loga,)x,y — (On/0 logX,)y,a] d logd, 
+ (On/0 loga,)x,y d log(a/)» 
Since 
(On/O logan)x,y — (On/0 logXx)y.a = 
[On/0 log(a/d)n)y.g — (On/O log(a/d)n)¥,a.......(61) 


symmetry suggests that this term be zero, so that as a function 
of Y and a/@, the differential of nm becomes 


dn = [(On/0 logXx)y.q + (On/0 log Vnsi)xa] d log Vnas 
+ (On/0 loga,)x,y d log(a/d)n 
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Justification for setting (61) equal to zero can be found by 
referring to a logX versus logY diagram. On this type of plot, 
lines of constant a and constant ¢ are parallel and of unit slope. 
A little manipulation will show that if log Y, — log Y,+41 changes 
from one tray to the next, then it makes no difference if this 
change is accomplished by a change in a or @, or both, the end 
result will be the same. 

If m had been expanded by means of the Taylor formula 
about loga,+:, and if Y,4; had been eliminated from the resulting 
expansion, there would be the term 


(On/0 logans1)x.y + (On/0 log Vait)xa = 
[On/0 log(Qn+1/On)) x6 — [On/0 log(an+1/hn)] xa... .(63) 


Following the same sort of reasoning that was employed in 
setting (61) equal to zero, this term can also be set equal to zero. 
There are then, the pair of equations 
[(On/O loga)x,y — (On/0 logX»)y,a] a=a, = 0 
[(On/0 loga)x,y + (On/0 log Yns1) x0] aani, = 

The difference in sign that appears here is due to the fact 
that A loga plays a dual role. Near X = Xn41, A loga behaves 
as — A logX, while near Y = Y,, it behaves as + A logY. 
Writing (64) and (65) for the same value of relative volatility 


(On/0 loga)x,y| = |(On/0 logXn)y,a| = |(On/0 log Vn+1)x,a|.. .(66) 


and 


(On/0 logXn)x.a = (On/0 log Yn41) x0 


Equation (67) simply states that total condensation of a part of 
the vapor stream and total vaporization of a part of the liquid 
stream yield the same results. 
Substituting (64) in (62), the total differential for 2 becomes 
dn = [(On/0 logXn)y.a + (On/0 log Yn41)x,a] d log Vng1 + 
(On/0 logXn)y,q d log(a/d)n 


When a/¢ is constant, Equation (4) becomes 
log( Ynui/V1) = — n log(a/g) 
and differentiation yields 
(On/0 log Yatilasg = — 1/log(a/)n........... 
Making the necessary substitutions in Equation (68) 
— log(a/@), dn = d log Yni1 + 3d log(a/@)n...........06. 


This is Equation (8) that was previously developed. 
In terms of liquid compositions, the companion equation is 


— log(Qnii/b,) dn = dlogX, — 3d log(an+1/Pn) 


If the equilibrium relationship had been written in terms of 
the equilibrium constant 


Yn ™ Kigktes. 
there would be the pair of equations: 
— log( Kuxn/yn+1) dn = d logynz1 + 3d log( Kyxn/¥n+1) 
— log( Kn+iXn/¥n41) dn = dlogx, — 3d log( Kn+iXn/¥n+1) 


The derivation of these equations parallels that of (71) since the 
geometry of the logy vs. logy diagram is the same as the log Y 
vs. logX diagram. 

The applicability of these equations is readily determined by 
evaluating the logarithmic composition terms, log( ¥2/Yn+1), 
log(Vn/¥n+1), ete., for adjacent trays and then dividing the 
arithmetic average of the two by the logarithmic average. The 
result gives the fraction of the actual trays accounted for by 
these equations at that point. This can be shown from Equa- 
tion (16) when log(a,/d,) = Qn is linear in log Y,41. Then, for 
one tray,“ 

4(On+ Dn) 
An = Ort FO) tog Qu1/Qn) 
(Qnt+1 — Qn) 


gives the calculated number of trays. When Qy+1/Q, lies between 
0.5 and 2, then Aw is less than 1.04. 


APPENDIX C: MATHEMATICAL CONVENTIONS 


A theoretical tray or plate is considered here as the difference 
operator E of the calculus of finite differences. E is such that 
E log Y, = logYn41. The subscript ” records the number of times 
the operator E occurs. Since E is related to the differential 
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Data 
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Comp 


operator d of the differential calculus by the Taylor series Basis: Assume (0n/0N), 
expansion; & = exp(d/dn); there is no difficulty with differ- . é 
ential theoretical trays. Accordingly, the differential tray is Equation (13) 

taken as the differential operator such that log Vasan = log ¥, + logio 2.06 


dlogY,. For a single theoretical tray it follows easily that lan 2.90/1.90; @, = 1.61 
10 Ms 


n+l 
(E — 1) logY, = d log Yn = log( Yn+1/Y,). Equation (27) 
o 1.61X, 


With these conventions, differential trays and differential numbers 1+1. 61, x, 
of theoretical trays are permissible. 


APPENDIX D: TRAY NUMBERING PROCEDURE 


lhe theoretical plate number, ”, can assume all values greater ; f 
than zero and its integral values need not coincide with a physical Equation (19) 
boundary. However, n itself represents an arbitrary boundary 4.6/0.3 
or reference plane. Thus, the quantity and properties of the ee == (5:75 7). 
reflux Sane 2 = and crossing a en plate boundary 1 + 24.6/0.3 _ L(x s + x4) 
bear the number of the plate as a subscript. Correspondingly, Oe 
the counter-flowing vane stream entering the plate Renles the (1.75/2.75) — a. 09/ 2.09) 24.6 + 0.3 
theoretical plate number increased by unity as a subscript. In "bie + x4) = 76.2 
effect, functional dependence upon 7 is indicated by a subscript ; : 
even though m is not a true index of demarkation. Equation (40) 
In this report, the counting of theoretical plates begins in (Lx), = 26.0/(20.6 — 1) = 1.3 
the enriching section at Xp = Yi; = Xo, and ends in the stripping (Lxs), = 9.0/(5.09 — 1) = 2.2 
section at Xz on the 45° diagonal of the McCabe-Thiele diagram ; 
in the equivalent binary system. The enriching section includes > . is ; 
the cuniaane and the stripping section includes the reboiler, Equation (42) & Feed Balance 
Equation (31) determines the theoretical plate requirements (Lxs)s = 11.0/(1 — 0.429/2.06) — 7 = 6.9 
for a separation from Xp to Xz. To determine the number of (Lx5); = 12/(1 — 0.206/2.06) — 10 = 3.: 
theoretical plates in the tower proper, it is necessary to deduct 
the equivalent number of plates for the reboiler and condenser. 
In the case of a total condenser and total reboiler the product 
compositions coincide with the terminal points of the tower ———_ 
operating lines and there is no correction. In the case of a partial C 
condenser, equilibrium reboiler, and thermosyphon reboiler, the . 
operating lines terminate before they reach the product composi- a 
tions at Xp and Xz. It is necessary to extrapolate the operating (Lx;) 
lines to evaluate the auxiliary equipment. In Equation (31), Mis 
this extrapolation was in effect accomplished by and evaluated Rateie’ 
by Equation (17). If correction is to be made to the total number a 
of plates obtained from Equation (31), it should be done by means 
of Equation (17). Customarily, the extrapolations are performed 
a ee the terminal flows at the top or Equation (18) (Enriching & Stripping Stepwise Application. ) 


The number of theoretical plates at the bottom of the enric hing . 1.75 ‘ 
section, e, and the number of plates at the top of the stripping — AN = | logi 24.6/0.3 blog io(2.06/3.12) 
section, f—1, are not necessarily equal unless the compositions atin dail 
of liquid and vapor at these points are equal. This is the case log 302.06 
in Equations (29) and (30) where saturated liquid and vapor 2.303 logic 
feeds are introduced at the intersection of the operating lines. logi03.12 

In the case of a partially vaporized feed in Equation (26), log (2.06/3.12) + 
the operating lines terminate before they intersect. A separation me 
effect results from the introduction of feed which is considered 0.4/16.7 
as a part of the distillation process and is included in the estima- log of - 175 
tion of tray requirements. In Equation (31), this composition oe 
difference is evaluated along extrapolations of the ope rating lines log101.86 
near their intersections. A number of different mechanisms can 2.303 log iol ——— 
be suggested for the coupling of the operating lines and the log 102.06 
introduction of partially vaporized feed. The simplest situation aa HR EP a 
occurs when he 6 ompositions of the liquid and vapor portions of logio( 1.86/2.06) 
the feed and the compositions at the intersection of the operating N = 10.91 (Enriching + Stripping) 
lines are the same and lie on the diagonal joining the equilibrium 
point (X»y, Yr) and the point Z on the 45° diagonal on the Equation (31) 
McCabe-Thiele diagram in the equivalent binary system. : 

n/N = 2.303 (2.90/1.90) logi2.90 = 1.62 


Reflux leaving the enriching section: 


+ $logio(1.86/2.06) 


n (including reb. & cond.) = 17.7; actual = 
APPENDIX E: SAMPLE CALCULATIONS Assume maximum flows are given by Equation (45): 


Data: Edmister, W. C., ‘‘Hydrocarbon Absorption and Frac- ‘ : 
tionation Process Design Methods”’, Tables 55, 58, 59. Reprinted Equation (45) 
from the Petrol. Eng. L.(x1 + Xu)e, max = 76. = 86.4:n <e 
Lil X1+-Xh mn, max > 76.2 2 + ( 22 +1) = 96.7; 
—<—— a = = — — ~ > f = 1 


| Dip | Bxe Fixe | ry : | Equation (19) 





86. 4(. X/1 +X)’, + 24.6 

V/1 + Y)"en. = = 

al dali tl a, 3644246403 °"~* 
96.1(X/1 + X)'n — 04 

V/1 + VY)" = ee fed 

i+ ' cntewees 


Intersection of the operating lines for maximum flows yields: 
Y"'41 = 1.63; X", = 1.05; 6", = 1.55 








The Canadian Journal of Chemical Engineering, June, 1961 137 





Equations (13) and (31) yield: 
n/N = 1.53; " = 16.7 


Inclusion in Equation (31) of the term neglected in the differ- 
entiation of Equation (14) results in: 


Equation (31a) 
N(On/ON )s log( On ON by 


n= mena 


+4] [(0n/ON) — 1] d log loga, 
(On/ON) . 
The correction term in Equation (31a) amounts to —0.32 in the 
first case and —0.24 in the second. 


Data: egg N. R., Pontinen, A. J., Ind. Eng. Chem., 
50, 730 (1958). — IV, 15th tes. Feed is vaporized on 
the feed tray: Xr = : Vp = = Jy. Cross plots of the data show 
that the intersection “f ‘the operating lines lies near the bottom 
of the feed tray. A check will be made of the number of trays 


Item z XB Xp x1 | Me xy W+1 
C, 0.200 | 0.009 | 0.853 | 0.910 | 0.903 | 0.233 | 0.367 
Cy | 0.370 | 0. 473 | 0.0165 | 0.042 | 0.0387 | 0.495 | 0.509 
a | | 2.00} 3.00 | 2.38 2.00 | 
Moles | 100 | 77.4 | 22.6 


(x; and y2 are compositions of liquid leaving and vapor entering 
the condenser. ) 

Vi41 2 Veer = 0.367/0.509 = 0.722 

X, = X;, = 0.233/0.495 = 0.470 

¢g; = 0.722/0.470 = 1.53 

(On/ON), = (logw2.00)/logi(2.00/1.53) = 2.60 


AN | 0.903 /0.0387 +H (2.38/3.00) 
— 2 ae 


> 303 | log 192.38 
£.9U9 10210) log 03.00 


~ Jogi(2.38/3.00) 


l eee + dogio(2.00/2.38) 
81" 6 903 /0.0387 cee 
:) 


log 192.00 
2.303 oun = 


+ 


log 192.3 


log(2.00/2.38) + 


, 0.009 /0.473 100 
Og 10! 0.722 /1OL 102. 


N = 10.75 (Partial Condenser + Enriching + Stripping) 
n/N = 2.303(2.60/1.60) logio2.60 = 1.55 


n (including reb. & cond.) = 16.7; actual = 16 


In the above example, heat of vaporization of the feed came from 
sensible heat in the feed. However, if the top tray and the feed 
tray act as heat transfer trays and heat effects are considerable, 
then stepwise application of Equation (17) may be necessary to 
avoid erratic results. Note that Equation (18) has been applied 
across the condenser itself because of the large change in relative 
volatility. 

The sum of the combined keys has a maximum of 152.0 mol./hr. 
in the liquid leaving the top tray, and a maximum of 154 mol./hr. 
in the liquid leaving the third tray beneath the feed tray. The 
operating lines for these flows are: 


: 152.0(X/1 + X)'» + 22.6(0.853) 
(F/t + Fyn ; 


152.0 + (0.853 + 0.0165 22.6 


(V/ + Vyfogy = BOM +X)" = 77.40.09) 
)nt1 = 154.0 — (0.009 + 0.473)77.4. eee 


Intersection of the operating lines yields, for the maximum flows: 


?”, 1.46; ”/N = 1.44; n” = 15.5 
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The calculated values of ” properly straddle the actual value. 
The correction term in Equation (3la) amounted to —0.2 in 


both cases. Some of the slopes of the operating lines in the 
equivalent binary system in arithmetic and logarithmic coordi- 
nates are shown below: 








AY/1 7 Y) \n+1 | A log Y; ntl 
n : Ea 
AC. X/1 + X)n A logX 
ee | Seca oe 
1 (cond.) 0.845 0.917 
2 (top) 0.886 | 0.923 
e 0.841 
Ss 0.848 0.909 


Let log(a/P)n = Qu — px where Q, is linear in log Y,41 and bn>0 
when the operating line lies to the left of the linear approximation 
and zero otherwise. Let log(a/®)n = Qn + gn where gn > 0 when 
the operating line lies to the Light. Then the operating line and 
the linear approximation will vield the same value of » when 


nd 1 ic. tee 
ie zg )aoern= |G, Grr) toe 9) 


Hence, | 
d0.| laud On| 
f" pes | | ee asda (84) 


In general, when pp/Qn? > gn/Q,? then the linear approximation 
will yield less than the actual number of trays. This is the 
situation which exists when the intersection of the operating 
lines lies to the left of the linear approximation in logarithmic 
coordinates. (It is understood in the above that dQ, changes 
sign at logY,41.) It is difficult to show that ” will always be 
less than actual when @,”’ is obtained from the maximum flows, 
although this appears to be the case. However, correction to n’’ 
for the linear approximation can be made. 


Data: Peiser, A. M., Chem. Eng., 67, No. 14, 129 (1960). 
Sufficient information is not available for a trav check. Design 
conditions will be analyzed. 





] — 
Comp. Xp Je | Se). Bp Xs | vy 
| 
vec lignan anil 
7 0.0592 | 0.0137 | 0.0002 | 0.1871 | 0.0379 | 0.0844 
8 0.1689 | 0.2243 | 0.0416 | 0.4351 | 0.3186 | 0.3877 
9 0.2759 | 0.1712 | 0.3963 | 0.0380 | 0.4440 | 0.3191 
10 0.2076 | 0.0706 | 0.2999 | 0.0003 0.1239 | 0.0797 
we 831.20 | 168.80 | 614.70 | 385.30 | 660.62 | 885.13 
salen ct at ee D site haaiiciaial 
Basis: | Components 8 and 9 will be taken as the keys. 
= 0.750 
Yr = 1.312 Y 41 = 1,223 Y, = 1.215 
Xp = 0.678 X, = 0.718 X;-1 = 0.70 
= 1.93 d. = 1.705 dy-1 = 1.735 


The intersection of the operating lines is estimated with Equa- 
tion (27): 


Vier = 1.185 
X, = 0.689 
$, = 1.72 = #(1.705 + 1.735) 


Using Equations (12), (31) and (44): 


(On/ON), = logy(1.93)/logi(1.93/1.72) = 5.70 
n/N = 2.303(5.70/4.70) logi(5.70) = 2.1 
(L/L); = 5.70/4.70 = 1.22 


Reverse fractionation in the stripping and enriching sections is 
avoided by preselection of the point of intersection of the operat- 
ing lines in the equivalent binary system. 


x. 7 .& 
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Factorial Design in the Study of Acid 


Leaching of Pegmatitic Uranium Ores’ 


D. G. FISHER?, R. G. McINTOSH:, R. L. EAGER® 
and A. B. VAN CLEAVE®* 


Factorial design of leaching tests done on pegma- 
titic uranium ores is an efficient method of studying 
the effects of changes in process variables, as the 
required information may be obtained with a mini- 
mum expenditure of effort. 


The results of a complete five factor, two level, 
factorially designed experiment in which an “acid- 
cure” technique was used are discussed. The data 
obtained in two, half replicate, five factor, two level, 
tests in which a percolation leach technique was used, 
establish base levels of variables for further economic 
study of uranium recovery from Higginson Lake 
pegmatites. 


large fraction of the known potential uranium reserves of 


the Canadian shield region of Northern Saskatchewan 
occurs in the form of pegmatitic deposits containing 0.10% or 
less of U30s. Mawdsley™ has described the geology of some 
of these deposits. The predominant minerals are feldspar and 
quartz but considerable amounts of biotite and smaller amounts 
of other minerals are also present. The uranium occurs mainly 
as small crystals of uraninite disseminated in the rock. The 
uranium values in these pegmatites show a tendency to concen- 
trate in the fines during comminution”?, which supports the 
suggestion that the uraninite is associated with the biotite or 
that it occurs along grain boundaries which are readily cleaved. 
Chemical and spectrographic analyses of Saskatchewan pegma- 
tites indicate that usually only trace amounts of carbonate, 
sulphide, phosphate, arsenate, molybdenum, cobalt, copper, zinc, 
nickel and chromium are present. Thorium in amounts compar- 
able to the uranium content is sometimes present. The low acid 
consuming character of the rock, as indicated by chemical 
analyses, has been confirmed by leaching tests which indicated 
that high uranium recoveries could be obtained by leaching the 
ore at 50 to 60% pulp density for periods of from 48 to 72 
hours with a solution containing 36 |b. of concentrated H.SO, 
per ton of crushed ore (<0.185 in. 23% minus 35 mesh) at 
ambient temperature and without the addition of an oxidizing 
agent), 
The results of leaching tests to be described in this paper 
were obtained using samples of a 15 ton lot of pegmatitic ore 
from the Anglo Barrington Mines property near Higginson 


1Manuscript received October 18, 1960; accepted March 18, 1961. 
2Carbide Chemicals Company, Montreal, Que. 

3Department of Chemistry and Chemical Engineering, University of Sas- 
katchewan, Saskatoon, Sask. 

Based on a paper presented at the C.1.C. Western Regional Conference, 
Regina, Sask. September 9, 1960. Contribution from the Department of 
ens —a = Chemical Engineering, University of Saskatchewan, Sas- 
atoon, Sas 
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Lake in Northern Saskatchewan. The average grade of the 
ore was 0.07% of U;Os. Preliminary leaching tests, in which 
several techniques were used with this ore have been reported. 

The purpose of this paper is to discuss the application of 
factorial design of experiments to the study of the variables 
affecting the recovery of uranium by acid leaching, using both 

“acid-cure” and percolation leach techniques. The problem is 
to evaluate the effects of a number of process variables; such as, 
time, temperature, amount of acid, acid concentration, oxidizing 
agent and particle size, on the recovery of uranium from a 
pegmatitic ore. The most economical investigation is made by 
arranging the experiment according to an ordered plan in w hich 
all factors are varied in a regular way. Provided that the plan 
has been correctly chosen, it is possible to determine not only 
the effect of each individual factor, but also the w ay in which 
each effect depends on the other factors, that is the interactions. 
Thus a more complete picture is possible than would be obtained 
by varying each of the factors one at a time while keeping others 
constant. Statistical analysis can be applied to the results of 
such designs and it is possible, if necessary, to have the design 
provide its own estimate of experimental error“. 

A factorial design enables one to obtain the required informa- 
tion with the required degree of precision and with a minimum 
expenditure of effort. If some previous knowledge is available 
regarding the effects of the factors and the interactions, it may 
be possibie to reduce the effort still further by using a fractional 
factorial design. The advantages of using factorial designs 
may be summed up as follows: 


1. When there are no interactions the factorial design gives 
the maximum effeciency in the estimation of the effects. 


te 


When interactions exist, their nature being unknown, a 
factorial design is necessary to avoid misleading con- 
clusions. 

3. In the factorial design the effect of a factor is estimated 
at several levels of the other factors, and conclusions hold 
over a wide range of conditions. 

The application of factorial design to “‘acid-cure” and perco- 
lation leach technique of recovering uranium illustrates the 
usefulness of the method in assessing the effects of process 
variables and in estimating the conditions that would result in 
maximum profit in a commercial operation. 


Experimental 
A. Acid-Cure Leaching 

In “acid-cure’”’ leaching the crushed ore is mixed with 
sufficient strong sulphuric acid solution to wet the pulp, which 
is then acid-cured for a period of time prior to washing to remove 
the pregnant solution. The procedure followed was: 
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TABLE 1 
Factor LEVELS For Actp CURE LEACHING EQUIPMENT 


400 g samples of crushed ore 





Factor Base | Unit —) Factor | | (+E ‘actor 
Level | Interval Level Level 

Temperature, °C. 30 10 20 40 
Amount of H2SO,, 

Ib/ton 27 9 18 36 
Volume of acid 

solution, ml. 60 | | 40 80 
Time, hr. 24 8 16 32 


<1.41 mm | <4.76 mm 
( —14 mesh )*|( —4mesh) 


Particle size 


“Sample II Table 2 >Sample I Table 2 

(1) 400 gm. samples of zrushed ore and the desired volume 
of sulphuric acid so.ution were added to sealers which 
were rolled at roon: temperature (23+3°C.) for 45 min. 
at 50 r.p.m. to ensure adequate mixing. 

(2) The samples were let stand in a constant temperature 
oven (t+0.5°C.) for specified times. 

(3) 200 ml. of wash solution (distilled water adjusted to 
pH 2 with H.SO,) was added and the sealer rolled for 
10 min. 

(+) The pregnant liguor was removed by suction filtration. 

(5) The residue was repulped in the orignal s ealer with 200 
ml. of wash solution by rolling for 10 min. 

(6) The residue was again filtered and repulped with another 
200 ml. of solution as before. 

(7) The residue was dried at 105°C. for 12 hours and sampled 
for uranium analysis. 

(8) The volumes of pregnant solution and the two wi ashings 
were measured and samples taken for uranium and R: Os 
analysis. The fluorophotometric method of analysis © 
was used for duplicate or triplicate analyses for uranium. 


A complete five factor, two level, factorial design was set 
up and the leaching tests performed in a random order. The 
factors and the levels chosen are summarized in Table 1. (Screen 
analyses for the samples used are shown in Table 2.) The levels 
were chosen to cover a range of convenient operating conditions 
and to — conditions which gave satisfactory uranium 
recoveries in previous experiments’ 39. The range between 
the upper a lower level of each factor was chosen so that 
the expected change in uranium recovery, produced by changing 
the level of one factor, would be approximately the same for 
each factor, thus minimizing the chance that the effect of 
changing one factor would over shadow the effects of other 
variables. ‘The range between the two levels of any factor was 
made large enough to ensure that expected differences in uranium 
recoveries would be large enough, when compared to the experi- 
mental error, to be tested effectively in the 32 tests. The factors 
and factor levels were chosen so that all treatments would be 
subject to the same experimental error. 


The per cent recoveries of U3Os in the leach liquor (including 
washings) for each of the 32 tests in the complete five factor 


design are shown in Table 3. The amounts of ignited R,O, 
that would be obtained from the neutralization of all of the 
leach liquor from 400 gm. samples of ore are also shown. These 
latter figures give an indication of the amount of iron, aluminum, 
silicon, etc. dissolved during the leaching process. It is evident 
that the greatest change in the weight of the ignited R.O; 
precipitates is produced by varying the amount of acid used 
and, that the amount of contaminants in the pregnant solution 
can be minimized by using the minimum amount of acid necessary 
for satisfactory uranium recoveries in a given time interval. 
Separate experiments indicated that between 8 and 20% fi of the 
acid soluble R»Os; fraction of the ore was dissolved in these 
leaching tests. 


An analysis of variance, by the method of individual com- 
parisons“ was done on the UsOsrecovery values coded according 
to the formula, ¢ = SinW% U;Os recovered. The statistical 
analysis showed that the main effects of the five factors were 
all significant at the 1% level. That is, it can be concluded, 
with only a 1 % chance of being wrong, that within the range 
studied, increasing the temperature, or increasing the amount 
of sulphuric acid or the volume of solution, or increasing the 
leaching time, or decreasing the average particle size, all result 
in a significant increase in the per cent of UsOxg extracted from 
the ore. 

The interactions between: (a) temperature and time of 
leaching; (b) amount of acid and particle size; (c) volume of 
solution and particle size; (d) amount of acid and volume of 
solution; were also significant at the 1% level. The interaction 
between temperature and volume of solution was found to be 
significant at the 5% level. The magnitudes of the various inter- 
actions are summarized in Tables 4 to 8 where the av erage per 
cent U;QOs extracted in the 8 tests* at the indicated levels of 
each pair of factors is shown. None of the other possible inter- 
actions were found to be significant. 


The results shown in Table 4 indicate that the time factor 
has a larger effect at 20°C. than at 40°C. The most probable 
reason for this interaction would be that the rate of leaching 
at 40°C. is so much greater than the rate at 20°C., that the 
extraction has proceeded practically to its full extent in 16 hours 
at the higher temperature level. If shorter leaching times has 
been selected this interaction might not have been significant. 


The amount of acid-particle size interaction shown by the 
data in Table 5 may be explained by assuming that 18 ib. of 
acid per ton of —14 mesh ore does not provide enough acid to 
fully extract the uranium, since a larger fraction of the acid 
would be consumed by reaction with the more finely divided 
rock. When 36 lb. of acid per ton of ore is used, the amount of 
uranium extracted is significantly greater with —14 mesh ore 
because more of the uranium mineral is in positions more 
accessable to acid attack in the time allowed. 


*Except in Table 7 where average U,O, recoveries for 4 tests are given. 


TABLE 2 
PARTICLE S1ZE DISTRIBUTION OF CRUSHED ORE SAMPLES USED IN ’’AcID-CURE”’ AND IN PERCOLATION LEACH TESTS 


Particle Size Interval 


mm, Standard Tyler Mesh. Sample I 
<4.76, >1.41 7 —4+14 | . 56.80 pcan 
<1.41, >0.23 | —144+65 | 31.92 | 
<0.23, >0.074 | —65 +200 | 7.58 | 
<0.074 | — 200 3.70 
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_Sample I 


| | 
16.43 | 6.8 11.3 
32 
| 


Percent 


Sample tit | Sample IV Sample V 


0.94 57.6 18.4 23.6 


75.31 34.4 | 64.7 | 58.1 
9.4 


— 
bo 
wn 
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oo 
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TABLE 3 


PERCENT RECOVERY OF U;03 AND AMOUNT OF IGNITED R2O3; IN 


400 gm. sample I (—4 mesh) 


LEACH LIQUORS FROM FACTORIALLY DESIGNED ’’ActD-CURE” LEACHING EXPERIMENTS 











aE, 40°C 
16 hr. 32 hr. 16 hr. 32 hr. 
H.SO, Vol. 7 i = ae: 
lb./ton ml, 
U,0 RO U;0; R»O; U;0s ROs U;0% RO; 
| : c g ; c g ’ c g ‘ Cc g 
36 40 84.3 1.96 87.2 2.50 87.0 2.86 87.3 2.92 
36 80 86.0 1.60 89.3 2.08 88.7 2.64 88.3 2.91 
18 40 84.0 1.89 86.0 1.54 83.5 1.65 86.8 1.58 
18 80 80.3 1.14 82.2 1.25 84.6 1.51 83.7 1.57 
400 gm. sample II (—14 mesh) 
20°C. 40°C. 
16 hr. 32 hr. 16 hr. 32 hr. 
H.SO, | Vol. a — a a nee - be es aaa al ce a 
Ib./ton | ml. 
U;0s RO; U,0 R»O U;0s ROs U;05 RO; 
% g % g % g Ye g 
36 40 87.9 2.30 91.0 2.86 92.0 3.17 91.2 2.96 
is i | _ 
36 80 89.9 a:i0 94.7 2.64 96.0 3.00 | 95.5 3.08 
18 | 40 83.6 Toe 86.2 1.66 85.2 1.71 85.1 1.60 
18 80 83.7 1.43 88.2 1.65 89.3 1.57 89.2 1.63 


The volume of solution-particle size interaction shown by 
the data in Table 6 may be explained as follows. With —4 
mesh ore, changing the volume of solution from 40 to 80 ml. 
resulted in a slight decrease in the amount of uranium extracted 
probably because the actual concentration of acid would be 


TABLE 4 


TEMPERATURE-TIME INTERACTION, 
“Actp-CurE” LEACHING 


Temperature 


20°C. 40°C. 
| 16 hr. 84.9 88.3 
Time | © of U;Os recovered 
| 32 hr. 88.2 | 88.8 
TABLE 6 


VOLUME OF SOLUTION-PARTICLE SIZE INTERACTION 
““Actp-CurE”’ LEACHING 


oy Si 
Volume of Solution 


| 


40 ml. 80 ml. | 
| —4mesh} 85.8 85.4 
Particle | | ©) U;0s of recovered 
size —14 mesh 87.8 | 90.8 
| 
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smaller with the larger volume of solution with resultant decrease 
in the rate of acid attack. With —14 mesh ore, however, the 
uranium extraction is considerably more efficient when the 
larger volume of solution is used. This could be due to incomplete 
wetting of the solid particles when 40 ml. of solution is used on 


TABLE 5 
AMOUNT OF ACID-PARTICLE SIZE INTERACTION 
‘**AcipD-CURE”’ LEACHING 


Amount of Acid 


| 18 Ib./ton} 36 Ib./ton! 


87.2 | 


—4mesh| 87.8 
Particle % of U3;0s recovered 
size —14 mesh 86.3 92.3 
TABLE 7 


‘TEMPERATURE-VOLUME OF SOLUTION INTERACTION 
““Actp-CuRE”’ LEACHING 


—4 mesh ore 


—14 mesh ore 


Temperature 


eG... 4G, 


20°C. 40°C. | 
| 87.2 88.4 | 86.1 | % Us0s 
| 90.1 92.5 


Volume of 40 ml. 85.4 


solution 


| | 
| 
80 ml. 84.4 86.4 | recovered 
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TABLE 8 


AMOUNT OF ACID-VOLUME OF SOLUTION INTERACTION 
“ActD-CURE” Le ACHING 


Amount of Acid 


18 lb./ton | 361b./ton 


40 ml. 5. 88.5 


| 
| q of U;08 


Volume of 
me recovered 


solution | 80 ml. 2 91.1 


—14 mesh ore where the surface area would be considerably 
greater than in the case of —4 mesh ore, where 40 ml. of solution 
were judged sufficient to effectively wet 400 gm. of ore. 


The data shown in Table 7 indicate that the major interaction 
between temperature and volume of solution is due to the results 
obtained with the — 14 mesh samples. The effect of temperature 
is probably more evident when the —14 mesh ore is treated 
with 80 ml. of solution than with 40 ml. of solution, because 
of more effective wetting of the particles by the larger volume 
of solution as was suggested i in explaining the volume of solution- 
particle size interaction above. 

The apparent interaction between amount of acid and volume 
of solution, shown by the data of Table 8, can also be explained 
as being due to the effect of volume of solution and the surface 
area of the ore particles as discussed above. 


In addition to providing data for testing the main effects of 


the five factors, this factorial design showed that satisfactory 
uranium recoveries are possible with “acid-cure’ ” leaching tech- 
niques on this type of ore. However, this technique may not 
be practical as a commercial process because of the large amount 
of energy required to mix the ore and the solvent. In contrast, 
percolation leaching can be adapted to the treatment of large 
tonnages of ore. Further experiments were designed to test the 
effects of possible process variables on the extraction of uranium 
by a percolation leach technique. 
B. Percolation Leaching 

The work involved in carrying out a complete factorial 
design comprising 2" tests when » factors are tested at two 
levels can frequently be lessened by the use of fractional factorial 
dseigns®.®, That is, it is possible to investigate the main effects 


and the more important interactions by using only a fraction of 


the number of tests required for the complete factorial experi- 
ment. This procedure is justified when it is known from previous 
tests that the effects of higher order interactions are negligible 
or when a sufficiently reliable estimate of experimental error 
has been established. However, if it is assumed that interactions 
of all orders are real, it follows that in a half replicate factorial 
design each effect is confused with another effect. Confused 
effects are termed “‘aliases’’ and it is imperative to determine the 
aliases for any proposed fractional factorial design in order to 
avoid confusion of important effects. Since it appeared that 
higher order interactions between leaching variables could safely 
be neglected, a half replicate factorial design was set up to test 
the effects of time, temperature, concentration of acid, concen- 
tration of oxidizing agent and particle size on recovery of 
uranium when a percolation leach technique is used. 

All percolation leach tests were conducted in a cabinet which 
permitted both temperature and humidity control. The air in 
the cabinet was continuously circulated by means of a fan. The 
leaching apparatus is shown diagramatically in Figure 1. Four 
such units were located in the cabinet. Crushed ore samples 
weighing 200 gm. were placed in the 4-in. Buchner funnels and 
the acid solution was forced up through the ore bed by means 
of a finger pump, at a rate which was slow enough to prevent 
ore particles s from being carried off in the overflow as it returned 
to the mixing reservoir. Evaporation losses were further mini- 
mized by tape sealing an inverted beaker over the funnel collect- 
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Figure 1—Schematie flow diagram of Percolation Leach 
apparatus. 
Legend: 1. 21. beaker; 2. Buchner funnel containing 200 g. 
of ore; 3. glass funnel; 4. Leach liquor in 600 ml. beaker; 
5. stirrer; 6. circulating pump. 


ing the overflow and by taping a flat cover over the mixing 
reservoir. Tests begun by placing the desired amount of leach 
solution in the mixing reservoir, bringing the entire system to 
the desired temperature and then starting the circulation of the 
leach liquor. At the conclusion of a test, the circulating pump 
was stopped and as much as possible of the leach liquor was 
drained into the mixing reservoir before being removed for 
measurement of volume and uranium analysis. The residue was 
washed by circulating 100 ml. of wash liquor (water adjusted 
to pH 2 with H.SO;) through the system for 14 hr. The first 
wash liquid was drained from the system aa the process 
repeated with a second 100 ml. portion of wash liquid. The 
volumes of both washings were measured and retained for 
uranium analysis. The leached ore sample was washed a third 
time by removing it from the Buchner funnel and pulping it 
with 100 ml. of distilled water. The residue was collected by 
suction filtration, dried, and analysed for uranium. The third 
wash liquid was also collected and analysed for uranium. 


The factor levels chosen for one half replicate factorial 
design I are shown in Table 9. Although earlier tests had not 
indicated the necessity of adding an oxidizing agent to obtain 
high uranium recoveries from this type of ore, it was of interest 
to see whether the presence of NaC1O; would significantly 
affect uranium recoveries, hence, the factor levels 0 and 4 Ib./ton 
of NaC1QO3 were chosen. The levels of other factors are similar 
to those used in previous experiments. After the experimental 
design had been chosen the 16 individual tests were randomized 
by drawing numbers from a hat. 

Screen analyses for ore samples III and IV used are given 
in Table 2. It will be noted that Sample III corresponds closely 
in size range fractions to Sample I, but that Sample IV shows 
a larger fraction of coarse particles than did Sample II, hence, 
it is designated as —8 mesh rather than — 14 mesh. 

The per cent uranium recoveries obtained in the tests of 
this half replicate factorial design are summarized in Table 10. 
In general, the uranium recoveries obtained by this leaching 
technique are very satisfactory, as high as 96% and the pregnant 
solution contains practically no suspended solids. Uranium could 
be recovered directly from such solutions by ion exchange or 
solvent extraction without employing a preliminary clarification 
step. 
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TABLE 9 
Factor LEVELS FOR PE RCOLATION- LE ACH De SIGN z 


200 gm. oiaiaiten of Crniied On 
me olume of Leach Solution 400, ml. 


Factor Base Unit |(—) Factor|(+) Factor 
Level Interval Level Leve 1 
Temperature °C. | 35 10 25 45 
Amount of H2SO, | 
Ib./ton | 2 9 18 36 
Amount of NaClO; | a 2 0 4 
Ib./ton 
Time hr. , 6 18 30 
Particle size <2.38 mm | <4.76 mm 
( —8me »sh )*|(—4 mesh)? 


‘aeneaie IV Table 2 >Sample III Table 2 


The responses were coded and subjected to statistical 
analysis by the method of Yates‘. The analysis of variance 
indicated that three of the main effects, temperature, amount of 
acid and amount of oxidizing agent were significant at the 1% 
level. The effect of time was significant only at the 5% lev el 
and the particle size factor was not significant under the condi- 
tions of this experiment. Only one interaction, that between 
temperature and particle size was found to be significant and 
that at the 5% level. Increasing the temperature from 25°C. 
to 45°C., the amount of H.SO, from 18 to 36 lb./ton of ore, 
and the addition of 4 Ib./ton of NaC10Os all significantly increase 
the extraction of uranium. Despite the observation that time 
appeared as significant only at the 5% level, it is certain that 
time is a controlling f factor in the leaching mechanism as could 
be substantiated by ‘choosing lower levels for this variable. The 
analysis does indicate that uranium recoveries would likely not 
be significantly increased by employi ing leaching times beyond 
30 hours. The lack of significance in the particle size factor 
indicates that virtually all ‘of the uranium oxide in the crushed 
rock is accessible to acid attack even in the —4 mesh material. 


TABLE 10 


PERCENT RECOVERY OF U;O03 IN LEACH LIQUORS FROM ONE-HALF 
REPLICATE FACTORIALLY DESIGNED PERCOLATION-LEACH, 
DEsIGN I 


200 gm. sample III (—4 mesh) 
_400 ml. ac id solution 


Zo°G; 


Hi 2SO, Ib./ton LN he .| 18 hr. 30 hr. 
89.5 | 87.7 





200 gm. somal IV (8 mesh) 
400 ml. | ac id solution 


IoC. 


| 
H.SO, lb./ton | NaC10; Ib./ton [18 hr. | 30 hr. | 18 hr. 30 hr. 


36 | 0 | 86.4 | 


r= [saa 


36 4 |} — | 95.3 | 96.1 
18 | 0 — | 83.9 | 88.0 


4 
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TABLE 11 


TEMPERATURE-PARTICLE SIZE INTERACTION 
PERCOLATION N LEACH DE SIGN I 





Particle Size 


—8 mesh | —4 mesh 





87.2 
Temperature 


93.1 


Although the temperature—particle size interaction is con- 
fused with the three factor interaction, amount of acid—amount 
of oxidizing agent—time, the latter has been assumed to be 
negligible. A two-way table of mean per cent U;QOs recoveries 
calculated from the data in Table 10 is shown as Table 11. 
Increasing the temperature from 25°C. to 45°C. has a larger 
effect on the —8 mesh material than on the —4 mesh material. 
Thus, in spite of the absence of particle size as a main effect, 
it appears that higher uranium recoveries would be obtained 
with more finely crushed ore at 45°C. An economic balance 
would, of course, have to be made to justify the high tempera- 
ture—fine particle size leaching method in practice. The revenue 
from increased recoveries would have to more than offset the 
increased grinding and heating costs to justify using these 
conditions. 

In preforming the tests for percolation leach design I, it 
was observed that the volume of solution used (400 ml. per 
200 gm. sample) was considerably larger than necessary and as 
pointed out above, particle size variation in the range inv estigated 
did not show a significant effect. Consequently a second half 
replicate factorially designed percolation leach experiment with 
factor levels as indicated in Table 12 was carried out. The 
particle size distribution of the ore sample V used in this design 
is given in Table 2. Note that the particle size distribution for 
Sample \ V is very similar to that for Sample 1V. The experi- 
mental procedure employed was the same as outlined above for 
fractional factorial design I. The per cent uranium recoveries 
obtained in the tests of design II are shown in Table 13. Analysis 
of variance for the uncoded per cent U;Qs recoveries for this 
design showed the effects of temperature, amount of oxidizing 
agent, and time to be significant at the 1% level, and the effect 
of amount of acid to be significant at the 5% level. Neither the 
volume of acid solution or any of the first order interactions 
were significant at these levels. 

As expected, the significance of the time factor is more 
evident in design II than in design I, because of the choice of a 
lower base level. The principal effect of i increasing the temper- 
ature or of adding an oxidizing agent appears to be to increase 
the rate of solution of the uranium mineral, so that a given 
uranium recovery can be attained in a shorter time if elevated 
temperature and an oxidizing agent are used. The same uranium 


TABLE 12 
FACTOR Le VELS FOR PE RCOL: ATION LEACH De SIGN Th 


© 200 gm. samples of crushed ore 


Factor | Base | 


Unit | =) Factor|(+) ) Factor 
| Level 


| Interval | Level | Level 


Temperature °C. | 45 
Amount of acid 
Ib./ton 
Amount of NaC10; 
Ib. /ton 


Time hr. 
Volume of acid 
solution, ml. 














TABLE 13 


PERCENT RECOVERY OF U;03 IN LEACH LIQUORS FROM ONE-HALF 
REPLICATE FACTORIALLY DESIGNED PERCOLATION-LEACH 
DESIGN IT 


200 gm. ‘anemia V(-8 sini’ 
125 ml. acid solution 


aa. 45°C. 


H2SO, lb./ton NaC 10; Ib., ton | 6hr. | 12hr.| 6hr. | 12 hr. 


18 0 


| 71.6 oo Src | 92.5 
18 4 | — | 8.0 | 88.0 |} — 
36 0 — £625 | 8B | _ 
36 4 4.0) — | - ea 
Sen, mempleV(—Smech) 
200 ml. acid solution 
asec. | 45°C. 


HH SO, Ib./ton | 


NaC10; Ilb./ton | 6hr. | 12 hr.| 6 hr. | 12 hr. 


18 i 0 | 80.2 | 82.9 
18 4 80.0 | | 89.0 
36 0 | 79.0) — | | 91.5 
36 4 | — | 886 | 87.5 | 


recovery can be achieved at lower temperatures, and without 
the presence of an oxidizing agent, if a longer contact time can 
be used. Since the base level for volume of acid solution was 
set lower for design II than for design I, the actual concentration 
of acid was sufficiently high for effective leaching at both levels 
of volume of solution, so variation in this factor did not prove 
to be significant. In practice, therefore, the least volume of 
solution that would allow effective circulation of the acid solution 
through the ore bed should be used. Considering the results of 
fractional designs I and II, it is evident that uranium recoveries 
of the order of 95% can be obtained by percolation leaching of 
—8 mesh ore at 45°C, with as little as 18 lb. of H:SO, per ton 
of ore for 24 hours, with oxidizing agent present. Similar 
recoveries can be obtained, if up to 36 Ib. of H. SO, per ton of 
ore is used without the addition of an oxidizing agent, although 
a somewhat longer contact period might be required. The high 
uranium recoveries and the “cry stal clear” pregnant solutions 
obtained in the percolation leaching of Higginson Lake pegma- 
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titic ore makes this method attractive as a possible commercial 
method for use with this type of ore. 

Conditions which would result in maximum U;Qs recovery 
could be arrived at by using the results of designs I and II to 
select points along the path of steepest ascent at which further 
tests should be done. It is evident that by using higher temper- 
atures, longer time, and larger amounts of acid in the presence 
of an oxidizing agent, that higher recoveries of U;Os would be 
obtained. However, that point at which the increasing costs 
exceed the profit achieved by increasing the recovery would 
soon be reached. Therefore, the profit associated with the 
recovery of uranium by any leaching method should be chosen 
as the response to be maximized now that approximate base 
levels for the significant factors have been established. Cost 
estimates to be associated with each of the significant factors, 
temperature, time, amount of acid, amount of oxidizing agent 
and particle size would have to be estimated to establish the 
profit obtainable at a given recovery level. ‘The above reported 
factorially designed experiments, using estimated profit as the 
response ‘could then be used as a a starting point to determine the 
levels of the process variables that would be used in plant 
practice. The value of factorial design in minimizing the number 
of tests necessary to establish the desired operating conditions 
is evident. 
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INDUSTRIAL SECTION | 


Reactions in a Fluidized Coke Bed 
with Self-Resistive Heating’ 


H. S. JOHNSON? 


When electrodes are inserted in a fluidized bed 
of conductive coke particles, an electric current 
can be passed with sufficient power to raise the bed 
to a high temperature. This device becomes a chem- 
ical reactor when the fluidizing gas is appropriately 
chosen. The reaction conditions which obtain in such 
a reactor are sometimes of especial value in pro- 
ducing favorable yields in chemical reactions. 


The design and operating characteristics of labor- 
atory scale units are described. These are usually 
Vycor tubes fitted with rubber stoppers through which 
electrical and gas connections are made. Even with 
such simple apparatus, it is possible to operate at 
1500°C. 

Several chemical reactions have been studied in 
detail using reactors of this type. Among these are 
some which show commercial promise, such as the 
high temperature reaction of ammonia and hydro- 
carbons to form hydrogen cyanide. 


fluidized bed of conductive particles can be heated by the 

passage of electric current between electrodes w hich dip 
into the fluidized bed, and in fact, such a system is a useful heat 
source, or chemical reactor. The physical arrangement is shown 
in Figure 1, 

When used as a chemical reactor, the coke bed is fluidized 
by the reactant gas or mixture, which is thereby heated very 
rapidly to the bed temperature. The main application of such a 
reactor is to the following types of chemical reaction: 


(1) reactions with substantial energy requirements. This 
applies equally to sensible and latent heat, but most 
especially to the higher temperature levels — say above 
1800°F. 


(2) reactions which are favored by rapid heating. 
(3) reactions which are favored by high temperature, either 
through thermodynamic or kinetic effects. 

The types of chemical reactions which can be carried out in 
such reactors are sometimes limited by the nature of the bed 
material. So far, only carbon has given satisfactory results, 
and carbon becomes reactive with various gases, oxygen, water, 
CO, and others, at moderately elevated temperatures. Thus, 
when it is proposed to include such gases as reactants, the 
existence of the reactions with carbon has to be considered. 

Previous references to reactors of this type are quite rare, 
and appear to be confined to the patent literature. A German 


1Manuscript received November 11, 1960; accepted February 17,1961. 
2Research Group Leader, Shawinigan Chemicals Limited, Shawinigan, Que. 
Based on a paper presented at the C.1.C. Chemical Engineering Con- 
ference, Quebec City, November 7-9, 1960. 
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patent to F. Winkler of I. G. Farbenindustrie in 1928 is of 
particular interest, but, for whatever reason, the proposal did 
not lead to commercial utilization. 


Design Considerations Applied to Laboratory Reactor 

(1) Electrode Entry: [¢ became apparent quite early that 
side-entering electrodes would not be operable. The reason for 
this is that, at elevated temperatures, the segment of reactor 
wall separating the electrode entry points becomes conductive, 
and is soon seriously damaged. it seems to be a general rule 
that the electrode should not enter at a hot zone, especially at 
points where it contacts the wall and the fluidized bed simul- 
taneously. It can also be concluded that, in a 2-electrode system, 
one electrode could enter at side or bottom, but the other must not. 


(2) Radiation Shield: [t is advisable to shield the electrode 
entry points from direct radiation from the fluidized bed. A 
radiation shield in the form of a false roof can be used, but when 
operating at high temperatures, it must be assumed that any 
electrical insulator will become conductive. The only com- 
pletely safe insulator in this region is a gap left between solid 
elements. For this reason, on the laboratory scale, radiation 
shields are made in segments supported by the electrodes and 
thermocouple well. (Figure 2) 

(3) Fluidization: For safe electrical operation, the coke 
particles must be in motion in the vicinity of the electrodes. 
If a continuous path of stationary particles exists between 
electrodes, local overheating, resulting in very large currents, 
frequently occurs. Good dispersion of the ‘fluidizing gas is 
therefore an important aid to steady operation. On the labora- 
tory scale, a porous-disk diffuser is very useful. When this is 
not practicable however, multi-jet dispersers or bubble cap 
types can be used. Materials of construction are porous carbon, 
graphite, and Vycor. 

Thermal insulation of laboratory reactors is kept to a mini- 
mum so that the Vycor tube will be cooler than the bed, which 
may operate at temperatures far above the quoted service 
temperature of Vycor. In most of our laboratory work, there- 
fore, the greater part of the input energy is lost to the air. With 
the combination of air- -cooling and the low thermal conductivity 
of Vycor, it is quite possible to operate a laboratory reactor at 
1500°C. using rubber stoppers for end closures. 

In commercial reactors, heat losses must be minimized and 
refractories must therefore be capable of withstanding tempera- 
tures equal to those in the fluid bed. Very important practical 
developments have been made in a pilot plant installation under 
the supervision of another research group at Shawinigan. These 
developments include features of design and material selection 
that make sustained high temperature operation feasible, and also 
certain process dev elopments connected with hydrogen cyanide 
manufacture. 
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Figure 1—Fluidized coke reactor with self-resistive heating. 


Electrical Properties 

Because of the motion of the conductive particles, the electric 
current in the electrofluid reactor is not entirely steady. More- 
over, the current depends upon the degree of fluidization. When 
the particles are not in motion, heat is generated in the path 
between the electrodes and is not readily dissipated, so that the 
path becomes very hot. Under these circumstances, local 
temperatures can become high enough to vaporize carbon, 
leading to agglomeration of coke particles. 

Under most conditions, Ohm’s law is approximately obeyed. 
This is best shown by connecting an oscilloscope to show the 
current-voltage relationships. The strict Ohm’s law relationship 
is represented by a straight line when current is plotted against 
voltage. Any curvature, or loop formation is a measure of the 
deviation from Ohm’s law. Sketches of some of the observed 
patterns are shown in Figure 3. It must be emphasized that in 
the normal operation of the electrofluid reactor, Ohm’s law is 
obeyed, and that the cases shown in Figure 3 include undesirable 
deviations from normal practice. Some understanding of these 
deviations clarifies the preferred operation. 

As the A.C. voltage is increased from zero, it is quite 
common to see the patterns shown in Figure 3 appear in suc- 
cession. In particular at the higher range of voltage (about 200 
volts per inch) alternation between type (2) and type (3) is 
quite common. Type (3) is considered to be characteristic of 
arcing as can be shown by withdrawing one of the electrodes 
from the fluid bed so that an arc is visible above the surface. 
In this situation the pattern is of type (3). It may be noted 


OHM'S LAW INTERMEDIATE ARCING 





Figure 3—Types of oscilloscope patterns observed in elec- 
trofluid reactor operation. 
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Figure 2—Forms of radiation shield used in laboratory 
reactor. 


that the average current under arcing conditions can be lower 
than under non-arcing conditions, even though the voltage is 
higher. If the average current, as measured on an ammeter, is 
plotted against average voltage, the above sequence is repre- 
sented as shown in Figure 4. 

At point A, there can be a sudden decrease in current. 
This is thought to represent the onset of arcing. Because of 
the lower current, the temperature will drop, the are is ex- 
tinguished and operation returns to the high-current level, 
whereupon the cycle repeats. 


Thermal Efficiencies 


Most fluid bed reactors are heated by combustion of fuel, 
often in a separate chamber from the reaction zone. Heat 
transfer from burning fuel has two characteristics which make 
it less efficient at high temperatures (above 1800°F.). 

(1) The rate of heat transfer depends on the temperature 
difference between the flue gas and the fluid bed. At 
temperatures above 2200°F. this constitutes a serious limita- 
tion on production rate. 

(2) When heat is being transferred at a high temperature level, 
only a relatively small fraction of the total sensible heat in 
the flue gas can be transferred to the fluid bed. Recovery 
of the remainder, by low-temperature heat transfer, e.g., 
raising of steam, is required if reasonably efficient fuel 
utilization is to be realized. 

These limitations of rate and efficiency of heat transfer at 
high temperatures do not apply to the electro-fluid reactor since 
the energy is fed directly to the reactant system, and the 
efficiency of electrical heating is not affected by the tempera- 


CURRENT | 





Figure 4—Ammeter-voltmeter measurements in electrofluid 
reactor (schematic). 
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ture level at which it is transferred, except in the minor effect, 
conimon to all types of heaters, of increased shell losses. 


Applications 

(1) CS, from petroleum coke and hydrogen sulphide. 
Available thermodynamic data show that the reaction between 
H.S and carbon is favored at high tem peratures ». By extra- 
polaticn it can be estimated that at 1327°C., 59.2% of the H:S 
fed reacts to form CSe, 26.6% reacts to form the elements, and 
the remainder is unchanged. Since carbon is a reactant, and the 
electrodes are of graphite, successful operatic mn obviously depends 
on a very large difference in reactivity between graphite and 
the carbon in the fluid bed. Using the coke derived in the 
petro ileum industry by fluid coking of residues, and briefly 
referred to as fluid coke, this condition is satisfactorily fulfilled. 


In a 38 mm. O.D. laboratory reactor, operating at tempera- 
tures in the vicinity of 1500°C., about 70% of the H:S was 
converted (1.e., not recovered) with apparent efficiencies to CS2, 
in the range 90-100%. These were thought to be higher than 
the true values because of conversion of an contained in 
the petroleum coke, forming additional CS, 


Experiments of longer duration were made in a 5-inch I.D. 
reactor at 1400°C. These indicated that graphite would be a 
satisfactory construction material. Conversion of H.S was 
estimated at about 55%, with efficiencies to CS, again in the 
range 92-94%. The low yield of elemental sulphur, compared 
with the durendynanic prediction, may be noted. 

Feed rate of H.S was approximately 5-6 s.c.f.m./cu. ft. of 
petroleum coke, in both the 34 mm. and the 5-inch reactors, 
but this can be varied within limits which depend chiefly on 
coke particle size. Production rates correspon¢ ling to the above 
feed rate are 15-20 Ib. CS./cu. ft. coke/hr. 


Che coke used in this work was commercial petroleum coke 
from a refinery fluid coker. It was usually screened to remove 
particles coarser than 20 mesh, and calcined in an open dish at 
temperatures above 900°C. in order to impart electrical con- 
ductivity. Actual power requirement for a 38 mm. laboratory 
a x at 1500°C. is approximately | kw., or pf aoe 
mately 10 amp. at 100 v. For a 51 mm. O.D. reactor at 1500°C. 
the power requirement is approximately 2.6 kw., or 13 amp. 
at 200 volts. 


(2) CS, from petroleum coke and sulphur. Commercial syn- 
thesis of CS, is usually carried out using wood charcoal, 
which is sufficiently reactive to permit operation at 700-800°C. 
Even at such temperatures, heat transfer imposes limitations on 
the production rate of an externally-fired reactor. In the electro- 
fluid reactor, heat input is not a problem, and it was therefore 
of interest to determine at what temperature fluid petroleum 
coke would react at a reasonable rate. 

It is known that the yield in this reaction is not greatly 
affected by temperature 3) and that an equilibrium vield of 
85 could be expected at 1200°C. Although the heat of reaction 
is close to zero at 1000°C., the sensible heat required to dis- 
sociate sulphur vapor and heat it to reaction temperature is still 
sufhcient to make the reaction suitable for this type of reactor. 

Sulphur can be fed to the : ttom of the reactor as vapor at, 
or above, the boiling point; or the liqn ud can be fed through a 
wide graphite pipe which, at a point near the bottom of the 
reactor, becomes hot enough to vaporize the sulphur. In the 
latter method, nearly all of the thermal energy involved 1 
supplied by the fluid coke reactor. 


It was found that : 3 100°C. the reaction proceeded rapidly, 
with a conversion of 70% being obtained in the laboratory. On 
a larger scale of ann where longer contact times ‘would 
be possible, higher conversions can be expected. Laboratory 
production rates were approximately 20 |b. 
coke/hr. 

In this temperature range, silicon carbide components are 
rapidly attacked. Again, graphite was found to be the most 
useful material of construction for the reaction zone. 
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INDUSTRIAL SECTIO 


(3) HCN from ammonia and hydrocarbons. 4 nearly ideal 
reaction for the electrofluid reactor is the synthesis of the highly 
endothermic HCN from ammonia and hy drocarbons, or ammonia 
and coke. The heat of formation of HCN (gas) at 77°F. is 
approximately 60,300 B.t.u./Ib. mole. Moreover, the reaction 
is favored by rapid heating. Further, since HCN itself is not 
thermody namically stable with re pect to its elements, the HCN 
must be removed from the reaction zone as quickly as possible. 
When all of these factors are considered, the electrofluid reactor 
is seen to be particularly well-suited . this synthesis. 

By way of orientation, it may be pointed out that the 
dominant route to HCN from aime is the Andrussow 
process“ for partial combustion of a mixture of CH, and NH; 
over a Pr catalyst at about 1000°C. Exit gases contain about 
8% HCN, and a typical yield on CH, is 70%, on NH3 60% 
Other high temperature processes using noble metal catalvsts 
are known. 

The object of the laboratory work was to evaluate the non- 
catalytic reaction between ammonia and hydrocarbons in which 
the fluid coke merely supplies the heat required. 

\ Canadian patent’® contains examples of operation at 
1400-1600°C., using propane and methane as hydrocarbons. 
In these experiments, the atomic ratio of nitrogen to carbon 
was varied, but was generally not far from unity. Gas feed 
rates are related to the size of coke particles use d, and in these 

experiments were in the range 5-10 volumes of feed gas (at 
room temperature) per volume of coke per minute. When the 
hydrocarbon is propane, the concentration of HCN in the off-gas 
from the reactor 1s about 25% by volume, unreacted ammonia 
about 0.6% by volume. The yield of HCN on ammonia fed is 
85-90% 
(4) HCN from ammonia and coke. This reaction is equally 
well-suited to the electrofluid reactor. Laboratory trials have 
shown that, in the absence of specific catalysts, ‘the yield of 
HCN is lower than when hydrocarbons are used as the source 
of carbon. 


Other Reactions 


The foregoing gives an idea of the methods used in studying 
certain reactions which have been of particular interest. In the 
five-year period since the development of the electrofluid reactor, 


numerous reactions have been studied. Some of these are merely 


listed: 

(1) COC 

(2) H:O-C 

(3) Chlorinations: 110.-Cl.-C 
$10.-Cl.-C 
\].03-Cl.-€ 


(+) Hydrocarbon cracking. 

(5) Metal oxides and sulphur, H.S, CS, to give metal sulphides. 

(6) Heating applications: The — 1 bed may be used as a 
preheater for chemical reactants, e.g., to preheat the recycle 
gas stream in iron reduction pr ocesses. 

(7) Py rolysis of acetic acid and acetone to give ketene. 
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